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ABSTRACT 


The  Pseudo  Wigner-Ville  Distribution  is  a  time-frequency  representation  of  an  input 
time  signal  and  is  ideally  suited  for  portraying  non-stationary  signals.  A  working  com¬ 
puter  program  is  presented  and  the  effect  of  preprocessing  and  postprocessing  data  ma¬ 
nipulations  is  shown.  The  program  has  been  developed  for  analyzing  data  for  use  in 
machinery  condition  monitoring  and  diagnostics  and  will  be  a  valuable  asset  for  ana¬ 
lyzing  transient  machinery.  A  practical  example  showing  pump  speed  variations  with 
time  is  also  presented.  Due  to  the  fact  that  the  Pseudo  Wigner-Ville  Distribution  can 
be  used  to  analyze  both  steady  state  and  transient  operations,  along  with  the  fact  that 
it  can  be  calculated  on  virtually  any  computer,  this  method  could  revolutionize  machin¬ 
ery  condition  monitoring  and  diagnostics. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  for  all  cases  of  interest.  While  even-  effort  has  been  made,  within 
the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic  er¬ 
rors,  they  cannot  be  considered  validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 


The  physical  condition  or  state  of  health  of  machineries  which  operate  in  short  du¬ 
ration  cycles  is  not  known  with  any  degree  of  accuracy.  Maintenance  on  these  ma¬ 
chineries  is  being  conducted  periodically  in  order  to  avoid  failures  and  prolong  the  useful 
operating  life  of  the  equipment.  These  machineries,  since  they  operate  for  only  short 
periods  of  time,  can  be  characterized  as  transient.  Additionally,  machineries  which  are 
not  operating  in  a  steady  state  condition  are  also  transient.  This  includes  rotating  ma¬ 
chinery  transitioning  from  one  speed  to  another. 

In  order  to  assess  the  physical  condition  of  machinery  without  complete  disassem¬ 
bly,  a  physical  measurement  of  its  vibrations  is  conducted  using  an  accelerometer.  Other 
sensors,  such  as  temperature  or  pressure  transducers,  could  also  be  used.  There  are 
other  methods,  including  motor  current  signature  analysis  on  electrically  driven  ma¬ 
chinery  and  wear  debris  analysis  which  could  be  used. [Ref.  1]  However,  vibrations  are 
used  predominantly  for  machinery  condition  monitoring.  The  vibrations  are  recorded 
in  the  time  domain. 

The  time  domain  representations  of  vibrations  may  be  decomposed  into  a  summa¬ 
tion  of  sine  waves  and  can  be  identified  in  the  frequency  domain  (see  Figure  1  on  page 
2  [Ref.  2]).  As  long  as  the  physical  characteristics  of  the  machinery  are  known,  these 
sine  waves  or  frequency  components  can  be  directly  attributed  to  physical  events  oc¬ 
curring  within  the  machinery.  For  example,  a  shaft  rotating  at  3600  RPM  (60  Hz)  with 
a  10  tooth  spur  gear  will  have  a  gear  mesh  frequency  of  600  Hz  (10x60  =  600). 
Therefore  a  frequency  component  of  600  Hz  may  be  attributed  to  the  gear.  As  the  speed 
of  the  shaft  changes  the  gear  mesh  frequency  will  also  change.  All  physical  components, 
including  bearings,  couplings,  impellers,  rotors,  etc.,  may  be  related  in  a  similar  manner 
to  frequencies  dependent  upon  the  shaft  rotation  speed.  Therefore,  for  transient  ma¬ 
chinery,  as  speeds  van.’  with  time  the  frequency  components  will  also  vary. 

There  is  a  need  for  a  method  to  represent  the  time  dependent  events  which  occur 
with  machinery  operating  in  transient  modes.  At  each  instant  in  time  as  the  speed  of  the 
machinery  changes  the  frequency  content  will  also  change.  The  Pseudo  Wigner-Ville 
Distribution  is  the  method  which  was  chosen  to  portray  these  time  dependent  changes. 
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Figure  1.  Time  and  Frequency  Domains 

The  Pseudo  Wigner-Ville  Distribution  is  a  three  dimensional  (time,  frequency,  and 
amplitude)  representation  of  an  input  signal  and  is  ideally  suited  for  portraying  transient 
phenomena.  The  Wiener  Distribution  has  been  used  in  the  areas  of  optics  [Refs.  3  ,  4, 
5]  and  speech  [Refs.  6,  7],  Wahl  and  Bolton  [Ref.  S]  are  using  it  to  identify  structure- 
borne  noise  components.  Flandrin  et.  al.  [Ref.  9]  recently  proposed  its  use  in  the  area 
of  machinery  condition  monitoring  and  diagnostics,  while  Forrester  [Ref.  10]  is  investi¬ 
gating  its  use  in  gear  fault  detection. 

The  Pseudo  Wigner-Ville  Distribution  can  be  used  to  portray  both  transient  non¬ 
stationary  phenomena  as  well  as  starionary  phenomena  and  therefore  can  be  used  for 
machinery  condition  monitoring  of  all  machinery.  This  includes  machinery  operating  in 
a  steady  state  condition.  Due  to  the  time  independent  nature  of  steady  state  operating 
conditions,  the  frequency  content  will  be  constant  for  all  time.  The  benefits  of  using  the 
Pseudo  Wigner-Ville  Distribution  and  monitoring  all  machinery  are  enormous.  Ma¬ 
chinery  which  has  never  before  been  monitored  now  can  be.  Additionally,  monitoring 
now  is  not  limited  to  just  a  given  steady  state  operating  condition,  all  speeds  can  be 
monitored,  the  effects  of  different  modes  of  vibration  can  be  investigated,  and  a  more 
thorough  evaluation  of  the  machinery  condition  can  be  obtained.  This  all  translates  into 
a  tremendous  economic  savings.  Here-to-fore  unmonitored  machinery  now  can  be 
monitored  and  machinery  monitoring  is  not  limited  to  just  steady  state  operating  con¬ 
ditions. 


The  WTGFL'N  computer  programs  presented  require  that  data  be  collected  in  the 
time  domain  and  be  digitized.  Once  it  is  digitized,  then  virtually  any  computer  can  take 
the  digitized  data  and  analyze  it.  The  only  hardware  required  is  a  transducer  with  a 
power  supply,  an  analog  to  digital  converter,  and  a  computer. 
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II.  THE  PSEUDO  WIGNER-VILLE  DISTRIBUTION 


A.  WIGNER  DISTRIBUTION  FUNCTION 

The  Wigner  Distribution  Function  (W'DF)  was  first  introduced  by  E.  Wigner  in  1932 
[Ref.  11].  Claasen  and  Mecklenbrauker,  in  a  three  part  series  of  papers,  developed 
mathematical  equations  for  the  WDF  in  continuous  time  [Ref  12]  and  discrete  time 
[Ref  13].  They  also  showed  relations  with  other  time-frequency  transformations  [Ref 
14],  For  the  continuous  time  case  using  two  different  complex  signals,  r(r)  and  s(/).  the 
cross-Wigner  Distribution  can  be  formed.  The  cross-Wigner  Distribution  is  defined  as: 

fOCe';coV(/  +  ^-)s*(/-f  )*  (1) 

— oe  * 

where: 

r  —  >•(/);  a  complex  time  signal 

s  =  sd):  a  complex  time  signal 

i  =  time 

cu  =  frequency 

*  =  complex  conjugate 

The  auto-Wigner  distribution  is  defined  as: 

IVDFS, ,(/,  co)  =  [' °°  e-j^s(,  +  )s'(t  -  ±  )dr  (2) 

—  DC  ^ 

Since  the  objective  of  this  research  is  to  accurately  portray  the  time  dependent  frequency 
characteristics  of  a  single  input  signal,  only  the  auto- Wigner  distribution  will  be  used. 
From  the  frequency  domain,  the  auto-Wigner  distribution  is  defined  as: 

V;DFS'S(oj,  0  =  -^rJ°C  e't'Siio  +  j-  )S*(cu  -  \  )d;  (3) 

where: 

5  =  S(u>)  ;  Fourier  Transform  of  s(i) 
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There  are  several  properties  of  the  WDF  which  are  important  to  note.  As  may  be 
seen  above,  the  WDF  in  the  frequency  domain  and  the  WDF  in  the  time  domain  are 
related  as  follows: 


U-DFSiS(co,r)=JrDFSiS(r,(o)  (4) 

A  time  shift  of  a  signal  is  a  time  shift  of  the  WDF: 

5(, _T)(/,  co)  =  ff'DFs  s(r  -  t,  co)  (5) 

A  frequency  shift  of  a  signal  is  a  frequency  shift  of  the  WDF: 

\VDFeia,s  e,a,s(t,  co)  =  1VDFSt s(t,  co -Cl)  (6) 

It  follows  that  a  time  and  frequency  shift  of  a  signal  is  both  a  time  and  frequency  shift 
of  the  WDF: 


]VDFe^),  e/%, -,)('•  ")  =  n'DFSt  s(t  -  T,  c O-Cl)  (7) 

Equation  (7)  is  extremely  important  when  we  consider  that  these  changes  in  time  and 
frequency  are  exactly  what  we  want  to  portray  for  our  use  in  characterizing  transient 
phenomena  for  machinery  condition  monitoring.  It  has  been  shown  that  the  WDF  can 
discriminate  the  frequency  content  of  a  signal  at  discrete  times. 

Integrating  the  WDF  over  time,  frequency,  and  both  time  and  frequency  provides 
signal  energy  information.  The  integral  of  the  WDF  over  frequency  at  a  specific  time 
yields  the  instantaneous  signal  power  at  that  time: 

J  irDFjJf,  co)dco  =  1 5(0 1 2  (S) 

The  integral  of  the  WDF  over  time  at  a  specific  frequency  yields  the  energy  density 
spectrum  of  a  signal  at  that  frequency: 

J""  1VDFS'S(t,co)dt=  |S(co)|2  (9) 
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The  integral  of  the  WDF  over  the  whole  plane,  both  time  and  frequency,  yields  the  total 
energy  in  the  signal: 


°°  f°°  lVDFSiS(t,  (o)dtd(o  =  |  |s|  |2  (10) 

-OO  —OO 

The  computer  program  listed  in  Appendix  A  does  not  accurately  maintain  the  energy 
information  in  the  original  input  signal.  In  the  preprocessing  program  the  input  time 
signal  is  windowed  and  amplified.  In  the  postprocessing  program  the  distribution  is 
averaged  over  both  time  and  frequency.  Additional  properties  of  the  WDF  are  listed  in 
references  12  and  15. 

The  discrete  time  auto-Wigner  distribution  as  developed  by  Claasen  and 
Mecklcnbrauker  [Ref.  13]  is  defined  as: 


U'DFS  s(t,  <o)  =  2  V  e~J2<°'s{t  +  t)s*(r  -  t)  (11) 

?=— DO 

The  time  and  frequency  shift  properties  of  the  continuous  time  WDF  described  by 
equations  (5).  (6),  and  (7)  remain  valid  for  the  discrete  time  WDF.  Yen  [Ref.  15]  defines 
the  discrete  time  auto-Wigner  distribution  for  a  sampled  time  signal  s(r),  for  0<  i<  T, 
as: 


T- 1 

1  I  .  4 -COT  , 

IF'DF,tJ(f.«)— ~ *'  +  *  ('"T)  (12) 

T=0 

Both  equations  (11)  and  (12)  are  similar.  In  either  case  the  WDF  is  basically  the  Fourier 
transform  of  an  auto  correlation  of  a  signal. 

B.  WIGNER-VILLE  DISTRIBUTION 

In  1948.  Ville  proposed  the  use  of  analytic  signals  in  time-  frequency  representations 
of  real  signals  [Ref.  16].  An  analytic  signal  is  a  complex  signal  which  contains  both  real 
and  imaginary  components.  The  advantage  of  using  an  analytic  signal  is  that  in  the  fre¬ 
quency  domain  the  amplitudes  of  negative  frequency  components  are  zero.  This  satisfies 
mathematical  completeness  of  the  problem  by  accounting  for  all  frequencies,  yet  does 
not  limit  our  practical  application  since  only  positive  frequency  components  have  a 
practical  interpretation. 


6 


An  analytic  signal  may  be  formed  from  a  real  signal  by  several  methods.  These 
methods  may  be  grouped  into  either  a  time  domain  formulation  or  a  frequency  domain 
formulation.  In  the  time  domain  a  typical  formulation  uses  the  Hilbert  transform  and 
may  be  expressed  as:  [Ref.  17] 

s{t)  =  sr(')+jH{sr (r)}  (13) 


where: 

s(/)  is  the  resulting  analytic  signal  which  is  complex 
s,(t)  =  real  component  of  a  signal 

=  imaginary  component  of  a  signal 

//{$,(')}  is  a  Hilbert  transform  and  is  defined  as: 


t  +  0 


t  =  0 


In  the  frequency  domain  a  typical  formulation  uses  the  Fast  Fourier  Transform  (FFT) 
and  is  expressed  as:  [Ref.  IS] 


s(t)  =  FFT-]lS(cj)-} 


(14) 


where: 


Si  w)  = 


SA(  cj)  cj  =  0 
2SR{a>)  W-1.---1 

0  otherwise 


s»  =  ffr[i,(/)] 


The  distribution  resulting  from  an  analytic  signal  being  processed  through  the  Wigner 
distribution  is  commonly  termed  a  Wigner- Ville  Distribution. 


C.  PSEUDO  W1GNER-VILLE  DISTRIBUTION 

The  Wigner- Ville  Distributions  of  most  signals  are  very  complicated  and  difficult  to 
interpret  since  the  input  signals  consist  of  many  components.  The  most  annoying 
characteristic  of  the  Wigner  and  Wigner- Ville  Distribution  is  the  presence  of  cross  terms. 
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Jones  and  Parks  [Ref.  19]  described  the  production  of  a  cross  term  for  the  sum  of  two 
signals.  The  Wigner  Distribution  of  the  sum  of  two  signals,  r{i)  +  s(i),  is  defined  as: 

WDFr+s<r+s{t,  to)  =  H'DFr  r(t,  o>)  +  2  Rc[)YDFr  s{t,  to)]  +  lYDFSiS(t,  to)  (15) 

The  annoying  cross  term  results  from  the  cross-Wigner  distribution  WDFri  and  is  lo¬ 
cated  midway  between  the  auto  terms.  As  the  input  signals  consist  of  a  summation  of 
greater  numbers  of  individual  components,  the  number  of  cross  terms  also  increases. 

There  are  several  approaches  for  the  removal,  or  deemphasis,  of  these  cross  terms. 
For  machinery  condition  monitoring  applications  the  presence  of  cross  terms  in  the 
Wigner- Ville  Distribution  is  not  disastrous  as  long  as  they  can  be  identified  as  such. 
However,  they  will  make  the  resulting  distributions  more  difficult  to  interpret.  Usually 
an  averaging  process  is  performed  in  order  to  produce  a  more  understandable  and  in¬ 
terpretable  representation  of  the  input  signal. 

There  are  two  methods  for  making  the  Wigner  Distribution  more  presentable. 
Claasen  and  Mecklenbrauker  describe  the  application  of  a  sliding  window  in  the  time 
domain  before  calculating  the  Wigner  Distribution  [Ref.  12].  The  resulting  distribution 
is  called  a  Pseudo  Wigner  Distribution.  A  second  option  is  to  smooth  the  Wigner  Dis¬ 
tribution  with  a  sliding  averaging  window  in  the  time-frequency  plane,  sometimes  re¬ 
ferred  to  as  a  Smoothed  Wigner  Distribution.  Since  the  effect  of  this  second  method  is 
essentially  the  same  as  the  first,  the  removal  of  undesired  components,  the  resulting 
distribution  from  either  method  will  be  called  a  Pseudo  Wigner  Distribution.  In  both 
cases  the  result  is  to  deemphasize  components  arising  from  calculations  and  to  empha¬ 
size  deterministic  components.  Obviously,  averaging  a  Wigner-Yille  Distribution  will 
result  in  a  Pseudo  Wigner-Yille  Distribution. 

There  has  been  a  fair  amount  of  research  as  to  the  optimum  smoothing  algorithm. 
References  12,  15,  20,  21,  22,  and  23  all  discuss  various  methods  available.  Due  to  the 
fact  that  the  resulting  Pseudo  Wigner-Yille  Distributions  were  to  be  used  to  identify  time 
varying  frequency  components  of  machinery  signatures,  a  sliding  exponential  window  in 
the  time-frequency  plane  is  used  in  this  work  for  the  averaging  process.  This  decision 
is  supported  by  Nuttall's  work  which  demonstrated  that  smoothing  must  be  applied  to 
both  time  and  frequency  [Ref.  23]. 
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III.  THE  W1GFUN  COMPUTER  PROGRAMS 


A.  BRIEF  DESCRIPTION 

The  computer  program  developed  in  this  research  is  included  as  Appendix  A.  The 
computer  code  begins  with  a  summary  of  the  programs,  subroutines,  and  symbols  used. 
The  computer  code  then  lists  the  W1GFUN  Fortran  programs,  is  followed  by  an 
alphnumeric  listing  of  the  subroutines,  and  is  concluded  with  a  listing  of  the  data  format 
conversion  program  referred  to  in  Appendix  B.  A  flow  chart  of  the  WIGFUN  programs 
follows: 


Figure  2.  Flow  Chart  of  WIGFUN  Computer  Programs 


The  program  is  written  in  Fortran  77  and  is  designed  to  be  interactive  so  few,  if  any, 
modifications  should  be  required  for  implementation  by  other  users.  The  programs  may 
not  have  been  validated  for  all  cases  of  interest.  The  'include'  statements  in  the  Fortran 
programs  cause  the  compiler  to  include  the  indicated  file  at  compilation.  The  included 
files  are  individual  subroutines  for  the  most  part,  making  the  size  of  the  individual  pro¬ 
grams  and  subroutines  more  manageable.  The  plotting  routines  all  use 
CA  -  D1SSPLA ®  software  [Ref.  24], 
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A  description  of  the  process  of  transferring  data  from  a  source  to  the  Digital  VAX 
computer  used  is  provided  in  Appendix  B.  The  only  requirement  for  the  data  to  be  used 
in  the  computer  program  is  that  it  is  in  the  format  read  in  subroutine 
DATA  1 N.  I XCLUDE. 

All  of  the  graphs  presented  in  this  thesis  resulted  from  analog  signals  obtained  in  the 
laboratory-  and  transferred  into  the  VAX  computer  using  the  procedure  provided  in 
Appendix  B. 

Figure  3  on  page  11  and  Figure  5  on  page  13  show  two  Pseudo  Wigner-Ville  Dis¬ 
tributions.  Figure  3  is  the  result  obtained  from  operating  on  two  sine  waves  generated 
by  function  generators  and  then  added  together.  It  can  be  clearly  seen  that  the  fre¬ 
quency  content  remains  constant  over  all  time,  as  expected  for  a  sine  wave.  Figure  4 
on  page  1 2  is  the  time  signal  input  for  Figure  3.  Figure  5  is  a  linear  chirp  which  was 
obtained  from  a  sine  wave  whose  frequency  was  linearly  varied  with  a  triangular  wave. 
This  figure  demonstrates  that  the  Pseudo  Wigner-Ville  Distribution  is  able  to  represent 
time  varying  frequencies.  Figure  6  on  page  14  is  the  time  signal  input  for  Figure  5. 
Figure  7  on  page  15  shows  different  viewing  angles  of  the  two  sine  waves  seen  in 
Figure  3.  The  top  left  view  shows  all  three  axes,  the  top  right  view  looks  down  the  fre¬ 
quency  axis,  the  bottom  left  view  looks  down  the  time  axis,  and  the  bottom  right  view 
looks  down  the  amplitude  axis  at  a  single  plane  located  at  one  third  of  the  maximum 
amplitude.  These  four  views  provide  a  good,  fast  presentation  of  the  distribution.  A 
more  detailed  contour  plot,  as  in  Figure  8  on  page  16,  gives  a  better  indication  of  the 
time  varying  characteristics  of  the  input  signal.  A  practical  machinery  monitoring  ex¬ 
ample  is  presented  in  chapter  4.  A  description  of  the  computer  program  listed  in  Ap¬ 
pendix  A  and  the  options  available  follows. 
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Figure  3.  Pseudo  Wigner-V'ille  Distribution  of  Tuo  Sine  Waves 


Amplitude 


Figure  4.  Time  Signal  for  Two  Sine  Waves 
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Pseudo  Wigner-Ville  Distribution 
Linear  Chirp 
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Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


e  5.  Pseudo  Wigner-Ville  Distribution  of  a  Linear  Chirp 
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Figure  8.  Detailed  Contour  Plot  of  Linear  Chirp 
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B.  LIMITING  DECISIONS 


As  can  be  seen  in  Figure  2  on  page  9,  the  first  Fortran  program  is  W1GFLX1. 
WIGFLX1  reads  in  the  raw  time  data,  preprocesses  it,  and  outputs  an  analytic  data  file. 
There  are  several  data  preprocessing  options  available  within  WIGFUNl  am.  each  will 
be  discussed.  However,  prior  to  preprocessing,  several  limiting  decisions  must  be  made. 

The  first  decision  is  the  number  of  data  points  to  use.  The  maximum  number  of 
data  points  which  may  be  processed  is  set  by  the  array  sizes  dimensioned  throughout  the 
computer  code.  As  printed  the  listing  is  set  for  512  data  points.  If  more  data  points  are 
required  then  larger  arrays  should  be  dimensioned.  Consequently,  more  computer  stor¬ 
age  space  is  required  and  the  run  time  is  longer.  It  was  found  that  for  initial  applications 
512  data  points  were  sufficient.  Due  to  the  use  of  a  Fast  Fourier  Transform  subroutine, 
the  number  of  data  points  must  be  a  power  of  2  (128,  256,  512,  1024,  etc.).  If  the  data 
file  to  be  used  does  not  have  as  many  data  points  as  were  chosen,  the  program  will  au¬ 
tomatically  pad  the  remaining  points  with  zeros.  The  program  will  interpret  this  as  an 
input  signal  value  of  zero. 

The  second  limiting  decision  which  must  be  made  is  what  time  step  size  ( A; )  to  use. 
Initially  the  program  reads  the  data  and  calculates  the  average  step  size.  From  this  av¬ 
erage  step  size  the  maximum  time,  maximum  frequency,  and  frequency  resolution  are 
calculated.  These  quantities  are  related  a<-  follows: 


nuime  —  dp  x  At 

(16) 

A/=  — - — — 

4  x  mlune 

(17) 

mfrcq  =  2  x  dp  x  A/ 

(IS) 

where: 

nuime  =  maximum  time 

mfieq  —  maximum  frequency 

dp  =  the  number  of  data  points 

At  =  time  step  or  time  resolution 

A/  =  frequency  step  or  frequency  resolution 

In  machinery  monitoring  applications  usually  two  frequency  spans  are  used,  a 
broadband  and  a  narrow  band  measurement.  These  displays  may  be  obtained  by  either 
changing  the  filters  before  digitizing  the  data  or  by  varying  the  time  step  size  in  the 
computer  program.  A  third  possibility,  which  has  not  been  implemented  but  could  be, 
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is  the  inclusion  of  a  digital  filter  in  the  computer  program.  By  selecting  the  smallest 
possible  At  the  largest  frequency  span  is  obtained.  The  computer  program  will  narrow 
the  frequency  span  and  focus  on  the  lower  frequencies  by  using  a  larger  At.  The  varying 
size  of  the  At  is  achieved  by  using  every  other  data  point  or  larger  multiples  of  data 
points.  Obviously  the  maximum  frequency  possible  is  determined  by  the  At  at 
digitization.  Sometimes  the  At  after  digitization  is  not  constant  and  the  program  allows 
for  this  by  calculating  an  average  At  and  using  this  average  value  as  the  actual  At  . 

C.  W1GFUN1 

As  mentioned  earlier,  W1GFUN1  is  a  time  domain  preprocessing  program.  It  reads 
in  the  raw  time  data  and  outputs  an  analytic  data  file.  Plotting  subroutines  are  included 
in  order  to  output  various  curves  throughout  the  process.  Before  the  analytic  signal  is 
calculated  there  are  several  options  which  may  be  exercised.  The  implications  and  ef¬ 
fects  of  each  are  discussed  below. 

1.  Effect  of  a  Mean  Value 

Figure  9  on  page  19  and  Figure  10  on  page  20  demonstrate  the  effect  of  a  mean 
value  in  the  time  domain.  Figure  9  is  a  sine  wave  whose  frequencies  have  been  slowly 
varied  and  which  has  a  mean  value.  As  can  be  seen,  the  mean  value  in  the  time  domain 
appears  as  a  DC  (0  Hz)  component  in  the  frequency  domain.  Figure  10  shows  the  same 
linear  chirp  which  has  had  the  mean  value  removed.  As  expected,  the  DC  component 
has  been  eliminated.  Figure  11  on  page  21  is  the  input  time  signal  for  the  linear  chirp 
shown  in  Figure  9  and  Figure  10. 
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Figure  10.  Linear  Chirp  with  Mean  Value  Removed  from  the  Time  Domain 
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2.  Effect  of  Amplification  of  Time  Signal  Amplitude 

A  subroutine  has  been  included  which  amplifies  the  amplitude  of  the  input  sig¬ 
nal  in  the  time  domain.  When  this  is  used  it  definitely  alters  the  energy  distribution 
representation.  The  only  other  effect  which  this  has  on  the  Pseudo  Wigner-Ville  Dis¬ 
tribution  is  that  it  serves  as  a  constant  multiplier,  vary  ing  the  amplitudes  of  the  distrib¬ 
ution.  Figure  12  on  page  23  and  Figure  13  on  page  24  depict  a  400  Hz  sine  wave  in 
noise  which  demonstrates  this  effect.  Figure  14  on  page  25  is  the  input  time  signal 
which  produced  Figure  12  and  Figure  13. 


Figure  12.  No  Time  Signal  Amplitude  Amplification 
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3.  Effect  of  Time  Signal  Windows 

A  subroutine  has  been  included  which  will  window  the  time  domain  signal.  The 
reason  for  windowing  the  signal  in  the  time  domain  is  to  bring  the  starting  and  finishing 
time  values  of  the  signal  to  zero.  There  are  two  windows  available,  a  modified  Hamming 
window  and  a  Hanning  window.  A  Generalized  Hamming  window  is  defined  as:  [Ref. 
25) 


window(t)  = 


a  +  (l  -ft)cos(-y-) 


0 


otherwise 


(19) 


A  Hanning  window  is  defined  as:  (Ref.  26] 


window(t)  — 


1  -  cos 


0  <t<T 


0 


otherwise 


(20) 


The  modified  Hamming  window  used  applied  a  cosine  weighting  function  to  the 
beginnning  and  ending  of  the  input  time  record  and  did  not  vary  the  amplitudes  inbe- 
tween.  The  equation  used  was: 

-^(l-cos-^y)  0  ^  r  <  0.1  T 


window(t)  = 


1 


o.ir</<o.9r 


1 

2 


^1  —  cos 


n{T  —  t)  \ 
0.17'  ) 


0.9  T<t<T 


0 


otherwise 


(21) 


A  modified  Hamming  window  is  the  preferable  window  since  it  alters  the  amplitude  of 
fewer  data  points  than  a  Hanning  window.  Figure  15  on  page  2S  shows  a  Pseudo 
Wigner-Yille  Distribution  of  a  2500  Hz  sine  wave  with  no  window.  Figure  16  on  page 
29  is  the  input  time  domain  signal.  Figure  17  on  page  30  shows  a  Pseudo  Wigner-Yille 
Distribution  of  a  2500  Hz  sine  wave  with  a  Hanning  window  applied  in  the  time  domain. 
Figure  IS  on  page  31  is  the  modified  time  domain  signal.  Figure  19  on  page  32  shows 
a  Pseudo  Wigner-Yille  Distribution  of  a  2500  Hz  sine  wave  with  the  modified  Hamming 
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window  applied  in  the  time  domain,  Figure  20  on  page  33  is  the  modified  time  domain 
signal.  In  Figure  15  it  can  be  seen  that  with  no  window  the  edges  of  the  graph  at  the 
start  and  finish  times  do  not  quite  reach  a  zero  value,  Figure  17  shows  that  the  Hanning 
window  drastically  changes  the  distribution,  and  Figure  19  shows  that  the  modified 
Hamming  window  provides  a  neater  distribution. 
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Figure  15.  Pseudo  Wigner-Ville  Distribution  with  No  Window 
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Figure  16.  Input  Time 
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Figure  19.  Pseudo  Wigner-Vilie  Distribution  with  Modified  Hamming  Window 
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4.  Effect  of  Making  a  Real  Signal  Analytic 

As  discussed  in  chapter  2,  an  analytic  signal  has  no  negative  frequency  compo¬ 
nents.  By  eliminating  these  negative  frequency  components  the  effects  of  aliasing  are 
greatly  reduced.  Figure  21  on  page  35  shows  the  aliasing  problem  which  results  when 
a  real  signal  is  not  made  analytic,  Figure  16  on  page  29  is  the  input  2500  Hz  sine  wave. 
Figure  22  on  page  36  shows  a  time  signal  which  has  been  made  analytic  and  has  had  a 
modified  Hamming  window  applied.  It  should  be  noted  that  the  window  was  applied 
before  the  signal  was  made  analytic,  thereby  saving  some  computation  time. 
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Figure  21.  Pseudo  Wigner  Distribution 
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Figure  22.  Analytic  Time  Signal 
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D.  WIGFUN1B 

WIGFl'Nlb  is  simply  a  program  for  plotting  the  analytic  time  signal.  It  is  set  ofT 
by  itself  in  Figure  2  on  page  9  since  it  is  not  required  to  be  run,  but  is  included  as  an 
added  convenience. 

E.  WIGFUN2 

W1GFUN2  reads  in  the  analytic  time  signal  from  a  file  and  calculates  the  Wigner- 
Ville  Distribution.  The  Wigner-Ville  Distribution  is  set  equal  to  the  real  pan  resulting 
from  the  complex  Fourier  transform  of  the  calculated  auto  correlation.  The  algorithm 
used  is  based  on  one  written  by  Wahl  and  Bolton  [Ref.  S]  and  can  be  expressed  as: 

v 

lVDFSJ{iAwjAt)  =  ^Re{FFT[2At  CORR(i )]}  (22) 


where: 

Re  =  Real  part 

CORR{i)  =  The  complex  auto-correlation 

CORR[i)  =  s(J  +  i  —  1  )s'{j  —  /  +  1)  1  <  /  <:  X  and  j  <  i 

CORRii  1  =  0  1  <  /  <  .V  and  j  >  i 

CORR(2.\  -  /  +  2)  =  CORR\i)  A'  <  /  <  2.Y 

1.  Effect  of  Changing  the  Definition  of  the  Wigner  Distribution 

As  shown  in  chapter  2.  there  are  several  different  definitions  possible  for  re¬ 
presenting  the  Wigner  Distribution.  The  computer  program  calculates  the  auto  corre¬ 
lation  of  the  signal  in  subroutine  CORR. INCLUDE  and  then  takes  the  complex  FFT 
of  the  result.  In  subroutine  WIGH. INCLUDE  the  WDF  was  set  equal  to  the  real  part 
resulting  from  the  FFT.  Other  possibilities  include  setting  the  WDF  equal  to  the  imag¬ 
inary  part,  setting  the  WDF  equal  to  the  square  of  the  real  part,  or  setting  the  WDF 
equal  to  the  square  root  of  the  sum  of  the  squares  of  the  real  and  imaginary  parts.  Ac¬ 
curate  representations  of  known  signals  were  obtained  using  just  the  real  part  so  this 
was  the  algorithm  used.  Additional  research  could  be  conducted  in  order  to  investigate 
other  options. 

F.  WIGFUN3 

WIGFUN3  is  a  postprocessing  program  with  three  functions.  It  reads  in  the  WDT, 
reduces  the  output  to  a  desired  size,  and  applies  a  sliding  averaging  window  to  the 
time-frequency  plane. 
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1.  Effect  of  Reducing  the  Wigner-Ville  Distribution 

If  the  Wigner-Ville  Distribution  is  not  reduced  and  one  started  with  512  data 
points,  the  resulting  distributions  would  be  1024  by  512  points.  Graphically  this  is  just 
too  many  points  to  visually  deal  with.  The  finest  resolution  which  is  understandable  is 
256  by  12S  points.  The  computer  program  allows  for  the  following  three  size  reductions; 
64  by  32.  12S  by  64,  and  256  by  128  points.  (See  Figure  23  on  page  39,  Figure  24  on 
page  40.  and  Figure  25  on  page  41.) 
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Pseudo  Wigner-Ville  Distribution 
Linear  Chirp 


16-APR-1990  09:50:29.07 
512  data  points 
Reduced  to  64  x  32 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplilied  by  0.0 
Hamming  window  time 


Figure  23.  Pseudo  Wigner-Ville  Distribution  Reduced  to  64  x  32  Points 
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Pseudo  Wigner-Ville  Distribution 
Linear  Chirp 


16-APR-1990  11:16:55.56 
512  data  points 
Reduced  to  126  x  64 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


Figure  24.  Pseudo  Wigner-Ville  Distribution  Reduced  to  128  x  64  Points 
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Pseudo  Wigner-Ville  Distribution 
Linear  Chirp 


16-APR-1990  12:17:40.91 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


Figure  25.  Pseudo  Wigner-Yille  Distribution  Reduced  to  256  x  128  Points 


2.  Effect  of  Smoothing  the  Wigner-Ville  Distribution 

As  discussed  in  chapter  2,  smoothing  the  Wigner-Ville  Distribution  deempha- 
sizes  the  terms  which  arise  from  mathematical  calculations  and  emphasizes  the  charac¬ 
teristic  events.  The  smoothing,  or  averaging,  process  also  makes  the  distribution  much 
more  understandable  (see  Figure  26  on  page  43  and  Figure  27  on  page  44).  It  is  im¬ 
portant  to  note  that  the  smoothing  process  uses  all  of  the  data  points  which  were  cal¬ 
culated  for  the  Wigner-Ville  Distribution  but  the  size  reduction  routine  only  plots  the 
size  desired.  This  averaging  also  effects  the  energy  representation.  Figure  28  on  page 
45  shows  the  different  views  and  Figure  29  on  page  46  shows  the  detailed  contours  of 
the  reduced,  but  not  smoothed,  distribution.  The  sliding  averaging  window  is  based  on 
an  algorithm  used  by  Wahl  and  Bolton  (Ref.  8  ].  At  a  given  point  on  the  time-frequency 
plane  the  weighted  window,  hg,  is  expressed  as: 

i  _jL  jL. 

hg{i,j )  = - - - e  200  200  (23) 

400-*AuA/ 

where: 

-  10  <  i  <  10 

-  10<j<  10 

The  weighted  window  is  applied  at  each  of  the  reduced  data  points  which  are  being 
plotted. 
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16-APR-1990  12:38:27.17 
512  data  points 
Reduced  to  256  x  128 


Mean  value  removed 
Time  amplified  by  0  0 
Hamming  window  time 


Figure  26.  Linear  Chirp  with  No  Smoothing 
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Pseudo  Wigner-Ville  Distribution 
Linear  Chirp 


16-APR-1990  12:17:40.91 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


Figure  27.  Linear  Chirp  with  Smoothing 
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16-APR-1990  12:24:05.01 
512  data  points 
Reduced  to  256  x  128 


Frequency  (Hz) 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


Figure  28.  Different  Views  of  Linear  Chirp  with  No  Smoothing 
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16-APR-1990  12:40:27.31 

512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  0.0 

_  _ Hamming  window  time 


Figure  29.  Detailed  Contours  of  Linear  Chirp  with  No  Smoothing 
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G.  WIGFUN4A,  WIGFUN4B,  AND  W1GFUN4C 

The  WJGFUN4  programs  plot  the  resulting  distributions.  Note  that  the  time  axis 
is  changed  from  seconds  to  milliseconds.  \\TGFUN4a  plots  the  64  by  32  point  distrib¬ 
utions,  WIGFUN4b  plots  the  12S  by  64  point  distributions,  and  \\TGFUN4c  plots  the 
256  by  128  point  distributions. 

H.  OTHER  CONSIDERATIONS 

1.  Effect  of  Noise 

Figure  12  on  page  23  and  Figure  13  on  page  24  bring  up  an  important  topic, 
which  is,  what  is  the  effect  of  noise  on  the  Pseudo  Wigner-Ville  Distribution? 
Figure  30  on  page  4S  is  the  Pseudo  Wigner-Ville  Distribution  of  two  sine  waves,  100 
Hz  and  500  Hz,  added  together  with  minimal  noise  present.  Figure  31  on  page  49  is  the 
input  time  signal.  Figure  32  on  page  50  is  the  Pseudo  Wigner-Ville  Distribution  of  the 
same  two  sine  waves  in  a  little  noise,  Figure  33  on  page  51  is  the  input  time  signal. 
Figure  34  on  page  52  is  the  Pseudo  Wigner-Ville  Distribution  of  the  same  two  sine 
waves  in  more  noise.  Figure  35  on  page  53  is  the  input  time  signal.  By  comparing  these 
distributions  it  can  be  noted  that  as  the  level  of  noise  is  increased  the  amplitudes  of  the 
distributions  van  but  the  frequency  content  remains  consistent.  Each  of  the  input  time 
signal  figures  are  made  up  of  an  upper  plot  from  W1GFUN1  (512  data  points)  and  a 
lower  plot  from  a  Dynamic  Signal  Analyzer  (1024  data  points)  recorded  at  digitization. 
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Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 


11-APR-1990  13:45:01.97 
512  data  points 
Reduced  to  126  x  64 
Smoothed  10  x  10 


30.  Pseudo  Wigner-Ville  Distribution  with  Minimal  Noise 
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Pseudo  Wigner-Ville  Distribution 
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11-APR-1990  14:30:35.97 
512  data  points 
Reduced  to  128  x  64 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
Hamming  window  time 
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Amplitude 


Pseudo  Wigner-Ville  Distribution 
100  Hz  and  500  Hz  sin  waves  in  noise 


11-APR-1990  14:36:00.50 
512  data  points 
Reduced  to  12B  x  64 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0.0 
HammlnQ  window  time 


Figure  34.  Pseudo  Wigner-Ville  Distribution  with  More  Noise 
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Figure  35.  Input  Time  Signal  «it!i  More  Noise 
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2.  Distribution  Interpretation 

The  interpretation  of  the  Pseudo  Wigner-Ville  Distribution  can  be  very  confus¬ 
ing.  Before  an  error  is  suspected  or  if  doubt  is  raised  as  to  the  accuracy  of  the  programs, 
verify  the  input  time  domain  signal.  The  Pseudo  Wigner-Ville  Distribution  will  accu¬ 
rately  portray  the  input,  however,  the  distribution  which  results  may  not  be  expected. 
Figure  36  on  page  55  is  the  distribution  from  the  input  time  signal,  Figure  37  on  page 
56.  Figure  38  on  page  57  is  the  distribution  from  the  input  time  signal,  Figure  39  on 
page  5S.  From  looking  at  the  input  time  domain  signals,  the  same  Pseudo  Wigner-Ville 
Distribution  could  be  expected,  however,  they  are  not.  The  peaks  in  Figure  38  result 
from  the  extra  data  before  and  after  the  chirp  in  the  time  signal. 
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Amplitude 
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Raw  Time  Signal 
Linear  Chirp  -  A 


3-  Swept  Sine  Wave  Example 

Figure  40  on  page  60  through  Figure  45  on  page  65  were  produced  from  a  sine 
wave  whose  frequencies  were  varied  with  a  sine  wave.  They  demonstrate  the  effect  of 
varying  the  Ai  and  the  ability  to  portray  time  dependent  frequencies.  Note  the  presence 
of  the  pronounced  peak  values  located  where  the  sweep  of  the  frequencies  are  changed 
from  positive  to  negative  and  negative  to  positive. 


mplitude 


Figure  40.  Swept  Sine  Wave,  Frequency  Span  0-1024  Hz 
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Pseudo  WIgner-Vllle  Distribution 
Swept  Sine  Wave 
Contours  from  zmax/6  to  zmax 


14-APR-1990  15:15:01.92 

512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  0.0 

Smoothed  10  x  10 _ _ Hamming  window  time 


Figure  41.  Detailed  Contour  Plot  Swept  Sine  Wave,  Frequency  Span  0-1024  Hz 
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Amplitude 


62 


Amplitude 


Pseudo  Wigner-Ville  Distribution 
Swept  Sine  Wave 


14-APR-1990  17:28:28.80 

512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  0.0 

Smoothed  10  x  10  _ _ Hamming  window  time 


Figure  43.  Swept  Sine  Wave,  Frequency  Span  0-512  Hz 
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14-APR-1990  17:30:43.08 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  0  0 
Hamming  window  lime 


Figure  44.  Detailed  Contour  Plot  S^ept  Sine  Wave,  Frequency  Span  0-512  Hz 
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IV.  PRACTICAL  EXAMPLES  USING  PSEUDO  W1GNER-V1LLE 

DISTRIBUTION 


The  figures  in  this  chapter  represent  data  taken  from  an  electrically  driven  single- 
stage  centrifugal  pump.  The  pump  speed  measurements  are  a  result  of  data  input  from 
a  proximity  probe  measuring  the  proximity  to  a  gear  directly  coupled  to  the  shaft. 

Figure  46  on  page  67  through  Figure  51  on  page  72  show  the  pump  operating  at  a 
steady  state  condition  with  flow  rates  of  10  gallons  per  minute  and  100  gallons  per 
minute.  Note  that  since  the  flow  rate  is  constant,  the  speed  of  the  pump  is  constant  for 
all  time.  The  speed  of  the  pump  for  high  and  low  flow  rates  is  not  very  different,  but  it 
can  be  seen  in  the  detailed  contour  plots  that  for  10  GPM  the  speed  is  about  52  Hz  and 
for  100  GPM  the  speed  is  about  4S  Hz. 

Figure  52  through  Figure  63  show  the  pump  in  transient  operations.  Figure  52  on 
page  73  through  Figure  55  on  page  76  show  the  pump  speeding  up  and  then  steadying 
out  at  a  constant  speed.  Figure  56  on  page  77  through  Figure  59  on  page  80  show  the 
pump  slowing  down  from  a  constant  speed.  Figure  60  on  page  81  through  Figure  63 
on  page  84  show  the  pump  operating  at  a  constant  speed,  slowing  down  after  being 
turned  off,  quickly  speeding  up  when  turned  back  on,  and  steadying  out  again. 
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Pseudo  Wigner-Ville  Distribution 
Pump  Speed  -  Steady  State  10  GPM 


14-APR-1990  18:49:30.16 

512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  5.0 

Smoothed  10  x  10 _ Hamming  window  time 


Figure  46.  Pump  Speed  Steady  State  10  GPM 
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Figure  47.  Detailed  Contours  of  Steady  State  10  GPM 
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Figure  48.  Input  Time  Signal  Steady  State  10  GPM 
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Pseudo  Wigner-Ville  Distribution 
Pump  -  Steady  State  100  GPM 


14-APR-1990  16:58:07.94 

512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  5.0 

Smoothed  10  x  10 _ Hamming  window  time 


Figure  49.  Pump  Speed  Steady  State  100  GPM 
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Pseudo  Wlgner-Vllle  Distribution 
Pump  -  Steady  State  100  GPM 
Contours  from  zmax/6  to  zmax 
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14-APR-1990  17:00:07.08 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  5.0 
Hamming  window  time 


Figure  50.  Detailed  Contours  of  Steady  State  100  GPM 
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Figure  51.  Input  Time  Signal  Steady  State  I0U  GPM 
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Pseudo  Wigner-Ville  Distribution 
Pump  Speed  -  Transient 
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14-APR-1990  17:16:26.47 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  5.0 
Hamming  window  time 


Figure  52.  Pump  Transient  Speed  Up 


Pseudo  Wlaner-VHIe  Distribution 
Pump  speed  -  Transient 
Contours  from  zmax/6  to  zmax 
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512  data  points  Mean  value  removed 

Reduced  to  256  x  128  Time  amplified  by  5.0 

Smoothed  10  x  10 _ Hamming  window  time 


Figure  53.  Detailed  Contours  of  Speed  Up 
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Pseudo  Wlgner-Vllle  Distribution 
Pump  Speed  -  Transient 


Mean  value  removed 
Time  amplified  by  5.0 
Hamming  window  time 


Figure  54.  Multiple  Views  of  Speed  t'p 
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512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 
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Pseudo  Wigner-Ville  Distribution 
Pump  Speed  -  Transient 


14-APR-1990  18:42:09.83 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  5  0 
Hamming  window  time 


Figure  56.  Pump  Transient  Coast  Don n 
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Figure  57.  Detailed  Contours  of  Coast  Down 
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Pseudo  Wlgner-VIlle  Distribution 
Pump  Speed  -  Transient 
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Figure  58.  Multiple  Views  of  Coast  Down 
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Pseudo  Wigner-Ville  Distribution 
Pump  Speed  -  Transient 


14-APR-1990  17:10:26.30 
512  data  points 
Reduced  to  256  x  128 
Smoothed  10  x  10 


Mean  value  removed 
Time  amplified  by  5.0 
Hamming  window  time 


Figure  60.  Pump  Transient  Coast  Down  and  Speed  Up 


Figure  61.  Detailed  Contours  or  Coast  Down  and  Speed  Up 
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Figure  62.  Multiple  Views  of  Coast  Down  and  Speed  Up 
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V.  CONCLUSIONS 


The  Pseudo  Wigner-Ville  Distribution  has  been  applied  to  analyzing  data  for  use  in 
machinery  condition  monitoring  and  diagnostics.  From  the  research  conducted  it  will 
be  a  valuable  asset  for  analyzing  transient  machinery.  The  following  conclusions  can 
be  drawn: 

•  The  Pseudo  Wigner-Ville  Distribution  is  ideally  suited  for  portraying  non¬ 
stationary  time  signals. 

•  Time  varying  frequency  components,  as  well  as  stationary  frequency  components, 
are  evident  in  a  single  output. 

•  The  use  of  an  analytic  signal  in  calculating  the  Pseudo  Wigner-Ville  Distribution 
eliminates  aliased  frequency  components. 

•  The  postprocessing  smoothing  process  deemphasizes  cross  terms  and  provides  a 
distribution  directly  related  to  the  input  signal. 

•  The  amplitude  of  the  Pseudo  Wigner-Ville  Distribution  varies  with  noise  and  with 
changing  frequencies  over  time. 

•  A  mean  value  of  an  input  time  domain  signal  is  a  DC  (0  Hz)  component  in  the 
Pseudo  Wigner-Ville  Distribution. 

•  A  modified  Hamming  window  applied  to  an  input  time  domain  signal  starts  and 
ends  the  signal  at  zero  amplitude  and  minimizes  the  distortion  of  the  Pseudo 
Wigner-Ville  Distribution. 

•  A  time  domain  signal  amplitude  amplification,  or  gain,  may  be  necessary  to  mini¬ 
mize  roundoff  error  in  calculating  the  Pseudo  Wigner-Ville  Distribution. 
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VI.  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


The  computer  program  works  well  and  the  Pseudo  Wigner-Ville  Distribution  will 
be  a  valuable  asset  to  assist  in  machinery  condition  monitoring.  Improvements  to  the 
computer  program  could  include: 

•  Increasing  the  number  of  data  points  being  processed 

•  Including  a  digital  filter 

•  Optimizing  the  routines  used  in  order  to  increase  the  speed  of  calculations 

•  Making  the  Pseudo  Wigner-Ville  Distribution  an  accurate  signal  energy  represen¬ 
tation 

Areas  where  additional  research  is  needed  include: 

•  Investigating  the  effect  of  using  just  the  real  part  of  the  Wigner-Ville  Distribution 
versus  the  real  and  imaginary  parts  or  combinations  thereof 

•  Investigating  other  sliding  averaging  windows  in  the  smoothing  routine 

•  Investigating  the  cause  of  the  pronounced  peak  amplitudes  evident  in  the  swept 
sine  wave  examples 

Recommendations  for  application  to  machinery  condition  monitoring  include: 

•  Using  an  analog  filter  on  all  input  data  to  ensure  unnamed  frequencies  will  not 
alias  the  frequencies  of  interest 
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APPENDIX  A.  COMPUTER  CODE 


A.  SUMMARY 
c 

c  [ rossano.  thesis] symbols,  include 

c 

c  This  is  a  description  of  the  programs,  subroutines,  and  symbols 

c  used  to  calculate  the  Pseudo  Wigner-Ville  Distribution 

c 
c 

c  FORTRAN  PROGRAMS 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c  SUBROUTINES 


c 


c 

AMPLIFY 

c 

Amplifies  the  time  signal 

c 

ANALYTIC 

c 

Converts  a  real  signal  into  an  analytic  one 

c 

CORR 

c 

Calculates  the  correlation 

c 

DATAIN 

c 

Reads  the  time  signal  data  into 

an  array 

c 

DATAOUT 

c 

Writes  the  WDF  array  to  a  file 

c 

DATA0UT2 

c 

Writes  the  RWDF  and  PWDF  arrays 

to  files 

c 

DTCALC 

c 

Calculates  delta  t  from  an  input 

file  or  array 

c 

FFT 

c  Calculates  the  Fast  Fourier  Transform 


WIGFUN1 

Reads  in  the  real  time  signal,  manipulates  it,  and 
writes  the  reculting  analytic  signal  into  a  file 

WIGFUNlb 

Plots  the  analytic  time  signal 

WIGFUN2 

Reads  in  the  analytic  time  signal,  calculates  the 
Wigner  Distribution,  and  writes  it  to  a  file 

WIGFUN3 

Reads  in  the  raw  Wigner  Distribution,  manipulates 
it  by  reducing  and  smoothing,  and  writes  the  results 
to  files 

WIGFUN4a 

Plots  the  3  dimensional  distributions  which  were 
reduced  to  (64,32) 

WIGFUN4b 

Plots  the  3  dimensional  distributions  which  were 
reduced  to  (128,64) 

WIGFUN4c 

Plots  the  3  dimensional  distributions  which  were 
reduced  to  (256,128) 
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a 


C 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


HAMMING 

Applies  a  hamming  window  to  a  signal 

HANNING 

Applies  a  harming  window  to  a  signal 

MAXAMP 

Finds  the  maximum  amplitude  of  a  signal 

MEAN 

Calculates  the  mean  value  of  a  signal 

MEANR 

Removes  the  mean  value  from  a  signal 

M INAMP 

Finds  the  minimum  amplitude  of  a  signal 

PL0T2D 

Plots  one  2  dimensional  curve 

PL0T2D2 

Plots  two  2  dimensional  curves  on  one  plot 

PL0T3D32 

Plots  one  3  dimensional  graph  of  rwdf(64,32) 

PL0T3D64 

Plots  one  3  dimensional  graph  of  rwdf( 128,64) 

PL0T3D128 

Plots  one  3  dimensional  graph  of  rwdf(256,128) 
PLOT3DSPLIT32 

Plots  four  3  dimensional  graphs  of  rwdf(64,32) 
in  a  split  screen  format  each  with  a  different 
viewing  perspective 
PL0T3DSPLIT64 

Plots  four  3  dimensional  graphs  of  rwdf( 128,64) 
in  a  split  screen  format  each  with  a  different 
viewing  perspective 
PL0T3DSPLIT128 

Plots  four  3  dimensional  graphs  of  rwdf(256, 128) 
in  a  split  screen  format  each  with  a  different 
viewing  perspective 
PL0TC0N32 

Plots  several  levels  on  a  single  contour  plot  for 
rwdf(64,32) 

PL0TC0N64 

Plots  several  levels  on  a  single  contour  plot  for 
rwdf( 128,64) 

PL0TC0N128 

Plots  several  levels  on  a  single  contour  plot  for 
rwdf(256, 128) 

PSEUDO 

Smooths  the  WDF  along  both  the  time  and  frequency  axes 

RANGE 3 2 

Finds  the  maximum  and  minimum  amplitudes  of 
rwdf(64,32) 

RANGE 64 

Finds  the  maximum  and  minimum  amplitudes  of 
rwdf( 128,64) 

RANGE 128 

Finds  the  maximum  and  minimum  amplitudes  of 
rwdf(256, 128) 

REDUCE 

Reduces  the  data  in  WDF  to  RWDF 
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SIZE 


TIMESIG 

WIGH 


WINDOW 


Changes  array  sizes  for  time  plotting  routines 
Modifies  and  plots  the  time  signal 
Calculates  the  WDF 

Calls  the  available  windowing  functions 


SYMBOLS 


ain( j )=signal  amplitude  array 

amp=input  data  amplitude 

ampl=amplif ication  applied  to  a  signal 

amplif=amplification  information  for  plot 

anq=answer  to  analytic  question 

atvq=answer  to  various  questions 

avgdelt=average  dt  found  in  input  time  array 

c(j)=complex  correlation  array 

ci( j )=imaginary  correlation  data  array 

coef=coef f icient  vsed  in  mathematical  calculations 

const=constant  used  in  mathematical  calculations 

cr(j)=real  correlation  data  array 

datetime=date  and  time  obtained  from  the  computer 

delt( j)=array  of  actual  dt  values 

delta=a  time  increment  used  in  calculating  hamming  window 

df=frequency  step 

dimf=dimension  size  for  frequency 

dimt=dimension  size  for  time 

dp=number  of  data  points 

dp2=dp"2 

dpend=modification  of  dp  in  input  array  gathering  to  ensure 
that  dp  number  of  points  will  be  collected 
dpnum=number  of  data  point  information  for  plot 
dt.=time  step 

dum~complex  variable  used  in  mathematical  calculations 

dum2=complex  variable  used  in  mathematical  calculations 

dum3=complex  variable  used  in  mathematical  calculations 

hg(i, j)=smoothing  multiplier 

i=a  counter 

ii=a  counter 

inname= input  file  name 

inv=indicator  for  fft  or  inverse  fft 

iplot_val=answer  whether  graph  output  on  screen  or  hardcopy 

j=a  counter 

jj=a  counter 

k=a  counter 

kk=a  counter 

L=a  counter 

LL=a  counter 

label=instructions  for  renaming  stdOOOOl.  dat 
lfile=plotting  file  name 
meani=mean  information  for  plot 
meanv=mean  value  of  a  signal 
mesh=mesh  size  for  plotting 
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c  mfreq=max  frequency  (Hz) 

c  mm=a  variable  used  for  reducing  the  number  of  time  points  in 

c  pwdf  and  rwdf  plots 

c  mnq=answer  for  mean  question 

c  mt=number  of  points  smoothed,  time 

c  mt2=mt*2 

c  mtime=max  time  (sec)  NOTE:  converted  to  (msec)  for  pwdf  and 

c  rwdf  plots 

c  n=a  counter 

c  nf=number  of  points  smoothed,  frequency 

c  nf2=nf*2 

c  nn=a  variable  used  for  reducing  the  number  of  frequency  points 

c  in  pwdf  and  rwdf  plots 

c  outname=output  file  name 

c  pi=mathematical  quantity,  pi 

c  pi2=2*pi 

c  pwdf(i, j)=pseudo  wigner  distribution 

c  rcnum=variable  for  which  correlation  to  plot 

c  rdp=reduced  number  of  data  points 

c  rdp2=rdp*2 

c  rnq=answer  for  rename  stdOOOOl. dat  question 

c  rval=reduction  value,  indicates  number  of  data  points 

c  reduced  to 

c  rwdf(i, j)=reduced  wigner  distribution 

c  s(j)=complex  time  array 

c  si( j)=imaginary  time  data  array 

c  sizedp=maximum  size  available,  set  by  dimension  of  arrays 

c  sr(j)=real  time  data  array 

c  ss=step  size  used  to  modify  dt  thereby  changing  max  freq 

c  startime=a  holder  for  datetime 

c  status=variable  used  with  datetime 

c  sum=sum  of  values 

c  sumb=sum  of  values 

c  sval=sine  of  val 

c  t(j)=data  array 

c  tcode=time  window  code 

c  tim=input  data  time 

c  tin( j )=signal  time  array 

c  title=title  of  a  plot 

c  titlehold=holds  3d  plot  title  name  while  plotting  correlation 

c  val=value  used  for  mathematical  calculations 

c  wdf(i,j)=raw  wigner  distribution 

c  windw=value  of  window  function  at  a  point 

c  wndwcode=type  of  window  code 

c  x(j)=data  array 

c  xmax=raaximum  amplitude 

c  xmin=minimum  amplitude 

c  y(j)=data  array 

c  ymax=maximum  amplitude 

c  ymin=miniraum  amplitude 

c  zhold=holder  for  calculating  zmax  or  zmin 

c  zmax=maximum  amplitude 

c  zmin=minimum  amplitude 

c 

c  NOTE:  The  plotting  routines  use  CA-DISSPLA  Version  11. 0 
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c 

c 


written  by  Computer  Associates  International,  Inc. 


B.  FORTRAN  PROGRAMS 

c 

c  [ rossano.  thesis] wigfunl.  for 

c 

c  Graham  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  is  the  first  of  four  that  calculate  and  plot  the 

c  Wigner  Distribution  Function  and  variations  of  a  signal 

c 

c  WIGFUN1  reads  in  the  real  time  signal,  manipulates  it,  and  writes 

c  the  resulting  analytic  signal  into  a  file 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  This  program  uses  an  include  statement  to  facilitate  changing 

c  the  size  of  the  plotting  arrays, 

c  (see  [ rossano. thesis] size,  include) 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  [ rossano.  thesis] symbols,  include 

c 


complex  s(1025) 

real  ain( 1025) ,mtime,mfreq,avgdelt ,dt ,sr( 1025) ,si( 1025) , 
&tin( 1025) ,df ,ampl 
character*23  datetime, tcode 
character*25  inname,  If ile 
character*50  title, outname 
character*80  label 
integer*4  status 

integer  dp,dimt,dimf ,rnq,dp2, j ,mnq 
include  ' size,  include' 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 

c  Filename  assignments 
print* 

print*,'  Program  to  calculate  and  plot' 
print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN1 ' 
print* 

print*,'  Enter  signal  input  filename' 

read(5,904)  inname 

print* 

print*,'  Enter  plotting  output  filename' 
read(5,901)  outname 
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print* 


c  Calculations 
dimt=1025 

call  dtcalc  (dimt , 1 , inname .tin.avgdelt) 
print*,'  Average  Delta  t  = 
write(5,906)  avgdelt 
print* 

print*,'  Input  the  number  of  data  points  to  evaluate' 

print*,'  (128,  256,  512,  or  1024)' 

print* 

print*,'  This  size  depends  on  the  size  of  the  arrays  set  in:' 

print*, '  [ rossano. thesis] size. include' 

print* 

print*,'  If  you  change  it  you  will  need  to  recompile  and  lnkdss' 
100  read(5,*)  dp 

print* 

if  (dp.  ne.  sizedp)  then 
print*,'  sizedp=' 
print*, sizedp 
print*,'  Reenter  dp' 
go  to  100 
endif 

if  (dp. eq. 128)  then 
dimt=130 
dimf=260 
endif 

if  (dp. eq. 256)  then 
dimt=260 
dimf=520 
endif 

if  (dp.  eq.  512)  then 
dimt=520 
dimf=1040 
endif 

if  (dp. eq. 1024)  then 
dimt=1025 
dimf=2050 
endif 

call  datain  (d  mt, inname, avgdelt, dp, ain,dt) 

tcode='No  window  time' 

wndwcode=0 

print*,'  Input  signal  has  been  put  into  an  array' 
print* 

call  timesig  (dimt,ain,dp,dt,outname,tcode, 

&title,s,sr,si, mnq , amp 1 ) 
open  (unit=7,file='wigl.  log' ,status='new' ) 
open  (unit=8 , f ile=' wigl. dat ' ,status=' new' ) 
write(7,908)  dimt,dimf ,dp,dt 
write(7,905)  outname.t code, mnq 
write(7,909)  title, ampl 
do  200  j=l,dp 
write(8,*)  s(j) 

200  continue 

close  (unit=7) 
close  (unit=8) 
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print*,'  Do  you  want  to  rename  std00001.dat?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  rnq 
print* 

if  (rnq.  eq.  1)  then 

print*.  Enter  Laser  plot  file  name' 

read(5,904)  lfile 

LABEL(1:  20)=' RENAME  STD00001.DAT  ’ 

LABEL(21:45)=lfile 
CALL  LIB$SPAWN( LABEL) 
print* 
end  if 

print*,'  Analytic  data  is  in  wigl.dat' 
print*,'  Plotting  info  is  in  wigl.log' 
print* 

print*,'  Now  run  WIGFUN2' 
print* 

c  Format  statements 

901  format( a50) 

904  format(a25) 

905  format(2x,a50,2x,a23,2x,i8) 

906  format( f 16.  8) 

908  format(2x,i8,2x,i8,2x, i8,2x,el6.  8) 

909  format(2x,a50,2x,fl6.  8) 
end 

c  SUBROUTINES 

include  'AMPLIFY.  INCLUDE' 
include  'ANALYTIC.  INCLUDE' 
include  ’ DATAIN.  INCLUDE’ 
include  'DTCALC.  INCLUDE' 
include  'HAMMING. INCLUDE' 
include  'HANNING. INCLUDE' 
include  'MAXAMP.  INCLUDE' 
include  'MEAN.  INCLUDE’ 
include  'MEANR.  INCLUDE' 
include  'MINAMP. INCLUDE' 
include  ' PL0T2D.  INCLUDE ’ 
include  ' PL0T2D2.  INCLUDE' 
include  ' TIMES I G.  INCLUDE' 
include  ' WINDOW.  INCLUDE' 
c 

c  [ rossano.  thesis] wigfunlb. for 

c 

c  Graham  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  plots  the  analytic  time  signal 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 


c  subroutines  may  be  found  in  [ rossano. thesis] symbols. include 

c 


complex  s(550) 

real  intime, mfreq,dt,sr(550) ,si(550) ,df ,ampl 
character*23  datetime, tcode 
character*25  lfile 

character*50  title, outname , t it lehold 
character*80  label 
integer*4  status 

integer  dp , dimt , dimf , rnq , dp2 , s izedp , mnq 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 
sizedp=512 

c  Get  data  from  WIGFUN1 
print* 

print*,'  Program  to  plot  the  analytic  time  signal' 
print* 

print*,’  WIGFUNlb' 
print* 

open  (unit=7 , file='wigl.  log' ,status='old' ) 
open  (unit=8,file='wigl.  dat' ,status=' old' ) 
rewind  7 
rewind  8 

read(7,908)  dimt , dimf , dp, dt 
read(7,905)  outname , tcode , mnq 
read(7,909)  title, ampl 
do  200  j=l,dp 
read( 8 ,*)  s(j) 

200  continue 

close  (unit=7) 
close  (unit=8) 
print*,'  Data  is  loaded' 
print* 

c  Calculations 

if  (dp.  gt. sizedp)  then 

print*,'  max  number  of  points  allowed  is' 

print*, sizedp 

print*,'  dp= 

print*, dp 

print* 

print*,'  Check  your  dimensions' 
stop 
endif 
dp2=dp*2 
mtime=dp*dt 
df=l.  /(4. *mtime) 
mfreq=2. *dp*df 
do  100  j=l,dp 
sr( j)=real(s( j)) 
si( j )=aimag( s( j ) ) 
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100  continue 
print* 

print*, '  Do  you  want  to  plot  the  analytic  time  signal?’ 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  ptaq 
if  (ptaq. eq. 1)  then 
titlehold=title 
title=' Analytic  Time  Signal' 
call  plot2d2  (sr ,si,dt,dp,title,outname,mnq) 
title=titlehold 
endif 

print*,'  Do  you  want  to  rename  std00001.dat?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  rnq 
print* 

if  (rnq.  ecj.  1)  then 

print*,  Enter  Laser  plot  file  name' 

read(5,904)  lfile 

LABEL(1:  20)=' RENAME  STD00001.DAT  ’ 

LABEL( 21: 45)=lfile 
CALL  LIB$SPAWN( LABEL) 
endif 

print*,'  WIGFUNlb.  FOR  is  complete' 

c  Format  statements 

904  format(a25) 

905  format(2x,a50,2x,a23,2x,i8) 

908  format( 2x, i8 , 2x, i8 ,2x, i8 ,2x,el6.  8) 

909  format(2x,a50,2x,fl6. 8) 
end 

c  SUBROUTINES 

include  ' MAXAMP. INCLUDE ' 
include  'MINAMP.  INCLUDE' 
include  ' PL0T2D2.  INCLUDE' 
c 

c  [  rossano.  thesis] wigfun2. for 

c 

c  Grahar!  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  is  the  2nd  of  four  to  calculate  and  plot  the 

c  Wigner  Distribution  Function  and  variations  of  a  signal 

c 

c  WIGFUN2  reads  in  the  analytic  time  signal  from  a  file, 

c  calculates  the  Wigner  Distribution  Function,  and  writes  the 

c  array  to  a  file, 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  {  rossano.  thesis] symbols,  include 

c 
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complex  s(550) ,c( 1100) 

real  mtime ,mfreq,dt ,cr( 1100),ci(1100),df, 
&wdf( 1100,550) ,ampl 
character*23  datetime, tcode 
character*50  title.outname 
integer,v4  status 

integer  dp,dimt ,dimf ,dp2,sizedp,mnq 
common  /wdfc/  wdf 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 
sizedp=512 

c  Get  data  from  WIGFUN1 
print* 

print*,'  Program  to  calculate  and  plot' 
print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN2' 
print* 

open  (unit=7,file='wigl.  log' ,status=’old' ) 
open  (unit=8,file='wigl.  dat' ,status=' old' ) 
rewind  7 
rewind  8 

read(7,908)  dimt ,dimf ,dp,dt 
read(7,905)  outname,tcode,mnq 
read(7,909)  title, ampl 
do  200  j=l,dp 
read(8,*)  s(j) 

200  continue 

close  (unit=7) 
close  (unit=8) 
print*,'  Data  is  loaded' 
print* 

c  Calculations 

if  (dp. gt. sizedp)  then 

print*,'  max  number  of  points  allowed  is' 

print*, sizedp 

print*,'  dp=' 

print*, dp 

print* 

print*,'  Check  your  dimensions' 
stop 
end  if 
dp2=dp*2 
mtime=dp*dt 
df=l. /(4. *mtime) 
mfreq=2. *dp*df 

call  wigh  (dimt,dimf ,s,c,df ,dt,dp2,dp) 

call  dataout  (dp) 

print* 

print*,'  Now  run  WIGFUN3' 
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o  o 


print* 


c  Format  statements 
905  format(2x,a50,2x,a23,2x,i8) 

908  format(2x,i8,2x,i8,2x,i8,2x,el6.  8) 

909  format(2x,a50,2x,fl6. 8) 
end 


c  SUBROUTINES 

include  ' CORR. INCLUDE' 
include  'DATAOUT.  INCLUDE’ 
include  ’FFT. INCLUDE' 
include  'WIGH. INCLUDE' 
c 

c  [ rossano.  thesis] wigfun3.  for 

c 

c  Graham  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  is  the  3rd  of  four  to  calculate  and  plot  the 

c  Wigner  Distribution  Function  and  variations  of  a  signal 

c 

c  WIGFUN3  reads  in  the  raw  Wigner  Distribution,  manipulates 

c  it  by  reducing  and  smoothing,  and  writes  the  resulting 

c  arrays  to  files 

All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 
c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  [  rossano.  thesis] symbols,  include 

c 


real  mtime,mfreq,dt,df , 

&wdf( 1100, 550), rwdf(256, 128), ampl 
character*23  datetime, tcode 
character*25  inname 
character*50  title, outname 
ir.teger*4  status 

integer  dp , dimt , dimf , dp2 , mnq , nn , mm , s izedp , rva 1 
common  /wdfc/  wdf 

c  Start  time  of  program 

status  =  lib^date_time  (datetime) 
print* 

type  *,  datetime 

sizedp=512 

nn=l 

mm=l 

c  Get  data  from  WIGFUN2 
print* 

print*,'  Program  to  calculate  and  plot' 
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print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN3 ' 
print* 

open  (unit=7,file='wigl.  log' ,status=' old' ) 
rewind  7 

read(7,908)  dimt,dimf ,dp,dt 
read(7,905)  outname , tcode ,mnq 
read(7,909)  title, ampl 
close  (unit=7) 

print*,'  Wigl.log  is  loaded' 
print* 


c  Calculations 

if  (dp. gt. sizedp)  then 

print*,'  max  number  of  points  allowed  is' 

print*, sizedp 

print*,'  dp= 

print* , dp 

print* 

print*,'  Check  your  dimensions' 
stop 
endif 
dp2=dp*2 
mtime=dp*dt 
df=l. /(4. *mtime) 
mfreq=2.  *dp*df 

open  (unit=13 , f ile=' wdf.  unf ' ,status=' old'  , 
&forir.=  ' unformatted'  ) 
rewind  13 

print*,'  reading  wdf  from  wdf. unf' 
print* 

do  100  j=l,dp 
do  200  i=l,dp2 

read(13)  wdf(i,j) 

200  continue 
100  continue 

close  (unit=13) 

print*,'  wdf  array  is  loaded' 

print* 

call  reduce  (dimt,dimf ,dt,dp,nn,mm,rwdf ,rval) 
print* 

call  dataout2  (rwdf ,3,dp,nn,mm) 
print* 

call  pseudo  (dimt.dimf ,rwdf ,dp,dt,nn,mm) 
print* 

call  dataout2  (rwdf ,2, dp, nn, mm) 
print* 

if  (title. eq. 'Wigner  Distribution')  then 
title=' Pseudo  Wigner  Distribution' 
endif 

if  (title. eq. ' Wigner-Ville  Distribution')  then 
title=' Pseudo  Wigner-Ville  Distribution' 
endif 

open  (unit=7 , f ile=' wig3. log' ,status=' new' ) 
write(7,908)  dimt ,dimf ,dp,dt 
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write(7,905)  outname ,tcode ,mnq 
write(7,907)  title, rval.ampl 
close  (unit=7) 

print*,'  Updated  plotting  info  is  in  wig3.log' 
print* 

if  (rval.eq. 1)  print*, '  Now  run  WIGFUN4a' 
if  (rval.eq. 2)  print*,'  Now  run  WIGFUN4b' 
if  (rval.eq. 3)  print*,'  Now  run  WIGFUN4c' 
print* 


c  Format  statements 
905  format(2x,a50,2x,a23,2x,i8) 

907  format(2x,a50,2x,i8,2x,fl6. 8) 

908  format(2x,i8,2x,i8,2x,i8,2x,el6.  8) 

909  format(2x,a50,2x,fl6. 8) 
end 

c  SUBROUTINES 

include  'DATA0UT2.  INCLUDE' 
include  ' PSEUDO.  INCLUDE ' 
include  'REDUCE.  INCLUDE' 
c 

c  [ rossano. thesis] wigfun4a. for 

c 

c  Graham  V.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  is  the  4th  of  four  to  calculate  and  plot  the 

c  Vigner  Distribution  Function  and  variations  of  a  signal 

c 

c  WIGFUN4a  plots  the  distributions  which  were  reduced  to 

c  64  x  32 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  [ rossano. thesis] symbols,  include 

c 


real  mtime ,mfreq,dt ,df ,rwdf(64,32) ,ampl 
character*23  datetime, tcode 
character*25  inname, If ile 
character*50  title , outname , titlehold 
character*80  label 
integer*4  status 

integer  dp , at vq , dimt , dimf , rnq , dp2 , s izedp ,n , mnq , 
&i , j , rval , rdp , rdp2 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 
sizedp=512 
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c  Get  data  from  WIGFUN3 
print* 

print*,'  Program  to  calculate  and  plot' 
print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN4a  reduced  to  64  x  32' 
print* 

open  (unit=7 , f ile=' wig3. log' ,status=' old' ) 
rewind  7 

read(7,908)  dimt,dimf ,dp,dt 
read( 7,905)  outname,tcode,mnq 
read(7,907)  title, rval,ampl 
close  (unit=7) 

print*,'  Wig3.log  is  loaded' 
if  (rval.  eq.  1)  then 
rdp=32 
else 
print* 

print*,'  You  are  using  the  wrong  program' 
if  (rval. eq. 2)  print*,'  Run  WIGFUN4b 
if  (rval.eq.  3)  print*,'  Run  WIGFUN4c' 
stop 
endif 

if  (dp. gt.  sizedp)  then 
print* 

print*,'  max  number  of  points  allowed  is' 

print*, sizedp 

print*,'  dp= 

print*, dp 

print* 

print*,'  Check  your  dimensions' 
stop 
endif 

c  Calculations 
dp2=dp*2 
rdp2=rdp*2 
mtime=dp*dt 
df=l. /(4.  *mtime) 
mfreq=2. *dp*df 

c  Convert  mtime  to  msec  for  plotting 
mtime=mtime*1000. 
print* 

print*,’  Do  you  want  to  plot  the  reduced  wdf?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq.  eq.  1)  then 

open  (unit=3, f ile=' rwdf . out ' ,status=' old' ) 
rewind  3 
do  100  j=l,rdp 
read(3,928)  n 
read(3,929) 
do  200  i=l,rdp2 

read(3,926)  rwdf(i,j) 


100 


200 

continue 

read(3,929) 

read(3,929) 

100 

continue 

close  (unit=3) 
titlehold=title 

if  (title. eq. ’Pseudo  Wigner  Distribution')  then 
t it le=' Reduced  Wigner  Distribution’ 
end  if 

if  ( title. eq. ’ Pseudo  Wigner-Ville  Distribution’)  then 
title=’ Reduced  Wigner-Ville  Distribution’ 
end  if 
print* 

print*,’  rwdf  is  loaded’ 
print* 

print*,’  Do  you  want  to  plot  the  rwdf  split  screen  view?’ 
print*,’  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dsplit32  (rwdf ,outname,mt ime, 
&mf req , dp , tcode , t it le , mnq , ampl ) 
print* 

print*,’  Do  you  want  to  plot  the  rwdf  single  view?’ 
print*,’  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3d32  (rwdf ,outname, 

&mt ime , mf r eq , dp , tcode , t it  1 e , mnq , amp 1 ) 
print* 

print*,’  Do  you  want  to  plot  the  rwdf  contours?’ 
print*,’  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotcon32  (rwdf ,outname, 

&mt ime , mf  r eq , dp , tcode .title, mnq , amp 1 ) 
ticle=titlehold 
endif 
print* 

print*,’  Do  you  want  to  plot  the  pseudo  wdf?’ 
print*,’  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 
if  (atvq. eq. 1)  then 

open  (unit=3,file=’pwdf.  out’ ,status=’ old’ ) 
rewind  3 
do  300  j=l,rdp 
read(3,928)  n 
read(3 , 929) 
do  400  i=l,rdp2 

read(3,926)  rwdf(i,j) 


400 

continue 

read( 3 , 929 ) 

read(3,929) 

300 

continue 

close  (unit=3) 
print* 

print*,’  pwdf  is  loaded’ 

print* 

print* 

print*,’  Do  you  want  to  plot  the  pwdf  split  screen  view?’ 
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print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dsplit32  (rwdf ,outname,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  pwdf  single  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3d32  (rwdf ,outname, 

&mtime ,mfreq, dp, tcode, title, mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  pwdf  contours?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotcon32  (rwdf .outname, 

&mt ime , rnf req , dp , tcode , t it le , mnq , amp 1 ) 
endif 

print*,'  Do  you  want  to  rename  std00001.dat?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  rnq 
print* 

if  (rnq.  eq.  1)  then 

print*,  Enter  Laser  plot  file  name' 

read(5,904)  lfile 

LABEL(1:  20)= 'RENAME  STD00001.DAT  ' 

LABEL(21:  45)=lfile 
CALL  LIB$SPAWN( LABEL) 
endif 
print* 

print*,'  Program  is  complete' 
print* 

c  Format  statements 

904  format( a25) 

905  format(2x,a50,2x,a23,2x,i8) 

907  format(2x,a50,2x,i8,2x,fl6.  8) 

908  format(2x,i8,2x,i8,2x,i8,2x,el6.  8) 

909  format(2x,a50) 

926  format( 2x,el6. 8) 

928  format(2x, i6) 

929  format(2x) 
end 

c  SUBROUTINES 

include  'PL0T3D3 2.  INCLUDE' 
include  ’ PL0T3DSPLIT32.  INCLUDE' 
include  'PL0TC0N32. INCLUDE' 
include  ' RANGE 3 2 . I NC LUDE ' 
c 

c  [ rossano. thesis] wigfun4b.  for 

c 

c  Graham  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 
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This  program  is  the  4th  of  four  to  calculate  and  plot  the 
Wigner  Distribution  Function  and  variations  of  a  signal 


c 
c 
c 

c  WIGFUN4b  plots  the  distributions  which  were  reduced  to 

c  128  x  64 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  [ rossano.  thesis] symbols,  include 


real  mtime,mfreq,dt ,df ,rwdf( 128,64) ,ampl 
character*23  datetime, tcode 
character*25  inname, lfile 
character*50  title,outname,titlehold 
character*80  label 
integer*4  status 

integer  dp , atvq , dimt , dimf , rnq , dp2 , s izedp , n , mnq , 
&i ,  j , rval , rdp , rdp2 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 
sizedp=512 

c  Get  data  from  WIGFUN3 
print* 

print*,'  Program  to  calculate  and  plot' 
print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN4b  reduced  to  128  x  64' 
print* 

open  (unit=7 , f ile='wig3.  log' , status=' old' ) 
rewind  7 

read(7,908)  dimt, dimf , dp, dt 
read(7,905)  outname, tcode, mnq 
read(7,907)  title, rval ,ampl 
close  (unit=7) 

print*,'  Wig3.log  is  loaded' 
if  (rval.eq.  2)  then 
rdp=64 
else 
print* 

print*,'  You  are  using  the  wrong  program' 
if  (rval.eq. 1)  print*,'  Run  WIGFUN4ar 
if  (rval.eq. 3)  print*,'  Run  WIGFUN4c' 
stop 
endif 

if  (dp. gt. sizedp)  then 
print* 

print*,'  max  number  of  points  allowed  is' 
print*, sizedp 
print*,'  dp= 
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print*, dp 
print* 

print*,'  Check  your  dimensions' 
stop 
end  if 

c  Calculations 
dp2=dp*2 
rdp2=rdp*2 
mt ime=dp*dt 
df=l.  /(4.  *mtime) 
mf req=2.  *dp*df 


c  Convert  mtime  to  msec  for  plotting 
mtime=mtime*1000. 
print* 

print*,'  Do  you  want  to  plot  the  reduced  wdf?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq.  eq.  1)  then 

open  (unit=3,file='rwdf.  out’ ,status=’ old' ) 
rewind  3 
do  100  j=l,rdp 
read(3,928)  n 
read(3,929) 
do  200  i=l,rdp2 

read(3,926)  rwdf(i,j) 

200  continue 

read(3,929) 

read(3,929) 

100  continue 

close  (unit=3) 
titlehold=title 

if  ( title. eq. ' Pseudo  Wigner  Distribution')  then 
title='Reduced  Wigner  Distribution' 
endif 

if  (title. eq. ' Pseudo  Wigner-Ville  Distribution')  then 
title=' Reduced  Wigner-Ville  Distribution' 
endif 
print* 

print*,'  rwdf  is  loaded' 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  split  screen  view?’ 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dsplit64  (rwdf ,outname, mtime, 
&rafreq,dp,tcode,title,mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  single  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3d64  (rwdf ,outname, 
&ratime,mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  contours?' 
print*,'  (1  for  yes,  2  for  no)' 
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read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotcon64  (rwdf .outname, 
Scmtime,mfreq,dp,tcode,title,mnq,ampl) 
title=titlehold 
end  if 
print* 

print*,'  Do  you  want  to  plot  the  pseudo  wdf?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq. eq.  1)  then 

open  (unit=3 , f ile='pwdf.  out' , status=' old' ) 
rewind  3 
do  300  j=l,rdp 
read(3,928)  n 
read( 3,929) 
do  400  i=l,rdp2 

read(3,926)  rwdf(i,j) 

400  continue 

read(3,929) 

read(3,929) 

300  continue 

close  (unit=3) 
print* 

print*,'  pwdf  is  loaded' 

print* 

print* 

print*,'  Do  you  want  to  plot  the  pwdf  split  screen  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dsplit64  (rwdf ,outname,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  pwdf  single  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3d64  (rwdf ,outname, 

&mtime ,mfr eq, dp, tcode .title ,mnq, amp 1) 
print* 

print*,'  Do  you  want  to  plot  the  pwdf  contours?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotcon64  ( rwdf .outname, 

&mt ime , mf  req , dp , tcode .title, mnq , amp 1 ) 
endif 

print*,'  Do  you  want  to  rename  std00001.dat?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  rnq 
print* 

if  (rnq. e^. 1)  then 

print*.  Enter  Laser  plot  file  name' 

read(5,904)  lfile 

LABEL(1: 20)= 'RENAME  STD00001.DAT  ’ 

LABEL( 21:  45)=lf ile 
CALL  LIB$SPAWN( LABEL) 
endif 
print* 
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o  o 


print*,'  Program  is  complete' 
print* 


c  Format  statements 

904  format(a25) 

905  format(2x,a50,2x,a23,2x,i8) 

907  format(2x,a50,2x,i8,2x,fl6. 8) 

908  format(2x,i8,2x,i8,2x,i8,2x,el6.  8) 

909  format(2x,a50) 

926  format(2x,el6.  8) 

928  format(2x,i6) 

929  format( 2x) 
end 

c  SUBROUTINES 

include  ' PL0T3D64.  INCLUDE' 
include  ' PL0T3DSPLIT64.  INCLUDE* 
include  ' PL0TC0N64.  INCLUDE' 
include  'RANGE 64.  INCLUDE' 

[ rossano.  thesis] wigfun4c.  for 
c 

c  Graham  W.  Rossano 

c  Prof.  J.  Hamilton 

c  Prof.  Y.  S.  Shin 

c 

c  This  program  is  the  4th  of  four  tc  calculate  and  plot  the 

c  Wigner  Distribution  Function  and  variations  of  a  signal 

c 

c  WIGFUN4c  plots  the  distributions  which  were  reduced  to 

c  256  x  128 

c 

c  All  of  the  subroutines  have  been  broken  off  for  ease  of  modifying 

c  and  printing.  Upon  compiling,  they  are  all  included, 

c 

c  A  description  of  the  symbols  used  in  this  program  and  its 

c  subroutines  may  be  found  in  ( rossano.  thesis] symbols,  include 

c 


real  mtime,mfreq,dt ,df ,rwdf (256 ,128) ,ampl 
character*23  datetime, tcode 
character*25  inname, If ile 
character*50  title ,outname,titlehold 
character*80  label 
integer*4  status 

integer  dp , atvq , dimt , dimf , rnq , dp2 , s izedp , n , mnq , 
&i , j , rval , rdp , rdp2 

c  Start  time  of  program 

status  =  lib$date_time  (datetime) 
print* 

type  *,  datetime 
sizedp=512 

c  Get  data  from  WIGFUN3 
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print* 

print*,'  Program  to  calculate  and  plot' 
print*,'  Wigner  Distribution  Function' 
print* 

print*,'  WIGFUN4c  reduced  to  256  x  128' 
print* 

open  (unit=7 , file=' wig3.  log' ,status='old' ) 
rewind  7 

read(7,908)  dimt.dimf ,dp,dt 
read(7,905)  outname,tcode,mnq 
read(7,907)  title, rval ,ampl 
close  (unit=7) 

print*,'  Wig3.log  is  loaded' 
if  (rval.eq.  3)  then 
rdp=128 
else 
print* 

print*,'  You  are  using  the  wrong  program' 
if  (rval.eq.  1)  print*,'  Run  WIGFUN’4ar 
if  (rval.eq. 2)  print*,'  Run  VIGFUN4b' 
stop 
endif 

if  (dp. gt. sizedp)  then 
print* 

print*,'  max  number  of  points  allowed  is' 

print*, sizedp 

print*,'  dp= 

print*, dp 

print* 

print*,'  Check  your  dimensions' 
stop 
endif 

c  Calculations 
dp2=dp*2 
rdp2-rdp*2 
mtime=dp*dt 
df=l.  / ( 4. *mtime) 
mfreq=2. *dp*df 


c  Convert  mtime  to  msec  for  plotting 
mtime=mtime*1000. 
print* 

print*,'  Do  you  want  to  plot  the  reduced  wdf?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq. eq. 1)  then 

open  (unit=3,file=' rwdf.  out' .status*' old' ) 
rewind  3 
do  100  j=l,rdp 
read(3,928)  n 
read(3,929) 
do  200  i=l ,rdp2 

read(3,926)  rwdf(i,j) 

200  continue 

read(3,929) 
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read(3,929) 

100  continue 

close  (unit=3) 
titlehold=title 

if  (title,  eq.  ' Pseudo  Wigner  Distribution')  then 
t it le=' Reduced  Wigner  Distribution' 
end  if 

if  (title,  eq.  ' Pseudo  Wigner-Ville  Distribution')  then 
t it le=' Reduced  Wigner-Ville  Distribution' 
endif 
print* 

print*,'  rwdf  is  loaded' 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  split  screen  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dsplitl28  (rwdf ,outname,mtime, 
&mf req , dp , tcode , tit le ,mnq , ampl ) 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  single  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plot3dl28  (rwdf ,outname, 

&mt ime, mf r eq, dp, tcode .title, mnq,ampl) 
print* 

print*,'  Do  you  want  to  plot  the  rwdf  contours?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotconl28  (rwdf ,outname, 

&mt ime , mf  r eq , dp , tcode , t it le , mnq , amp 1 ) 
title=titlehold 
endif 
print* 

print*,'  Do  you  want  to  plot  the  pseudo  wdf?' 
print*,'  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 
if  (atvq. eq. 1)  then 

open  (unit=3,file='pwdf. out' ,status='old' ) 
rewind  3 
do  300  j=l,rdp 
read(3,928)  n 
read( 3,929) 
do  400  i=l,rdp2 

read(3,926)  rwdf(i,j) 

400  continue 

read( 3,929) 
read(3,929) 

300  continue 

close  (unit=3) 
print* 

print*,'  pwdf  is  loaded' 

print* 

print* 

print*,'  Do  you  want  to  plot  the  pwdf  split  screen  view?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
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if  (atvq. eq.  1)  call  plot3dsplitl28  (rwdf ,outname,mtime, 
&mf req, dp, tcode, title, mnq, amp 1) 
print* 

print*,’  Do  you  want  to  plot  the  pwdf  single  view?’ 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 

if  (atvq. eq.  1)  call  plot3dl28  (rwdf ,outname, 

&mt ime , rof r eq , dp , tcode , t it le , mnq , ampl ) 
print* 

print*,'  Do  you  want  to  plot  the  pwdf  contours?' 
print*,'  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  plotconl28  (rwdf ,outname, 

&mt ime , mf  r eq , dp , tcode , t it le , mnq , amp 1 ) 
end  if 

print*,'  Do  you  want  to  rename  std00001.dat?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  rnq 
print* 

if  (rnq. eq. 1)  then 

print*,  Enter  Laser  plot  file  name' 

read(5 ,904)  lfile 

LABEL(1:  20)= 'RENAME  STD00001.DAT  ' 

LABEL(21:  45)=lfile 
CALL  LI B$ SPAWN (LABEL) 
end  if 
print* 

print*,'  Program  is  complete' 
print* 

c  Format  statements 

904  format(a25) 

905  format(2x,a5C,2x,a23,2x,i8) 

907  format(2x,a50,2x,i8,2x,fl6. 8) 

908  format(2x,i8,2x,i8,2x,i8,2x,el6.  8) 

909  format(2x,a50) 

926  format(2x,el6. 8) 

928  format(2x,i6) 

929  format(2x) 
end 

c  SUBROUTINES 

include  'PL0T3D128.  INCLUDE’ 
include  ' PL0T3DSPLIT128.  INCLUDE' 
include  ’ PLOTCON128.  INCLUDE’ 
include  'RANGE 12 8.  INCLUDE' 

C.  ALPHANUMERIC  LISTING  OF  SUBROUTINES 

c 

c  [ rossano. thesis] amplify,  include 

c 

c  This  subroutine  amplifies  the  signal  ain 


subroutine  amplify  (dimt,ain,dp,ampl) 
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integer  dp,j,dimt 
real  ain(dimt) ,ampl 

print*,  '  What  value  would  you  like  to  amplify  the  signal  by?' 
print*,  '  (Do  not  forget  the  decimal  point)' 

read(5,901)  ampl 
do  100  j=l,dp 

ain( j )=ampl*ain( j ) 

100  continue 
return 

901  format( f 16. 8) 
end 
c 

c  [ rossano.  thesis] analytic. include 

c 

c  This  subroutine  converts  a  real  signal  into  an  analytic  one 


subroutine  analytic  (dimt,sr,s ,dp) 

integer  dimt,dp,i,j 
complex  s(dimt) 

real  sr(dimt) ,dt,pi,sum,val ,n,sumb,sval 

pi=4.  *atan(  1.  ) 

do  100  i=l,dp 
sum=0. 

do  200  j=l,dp 
sumb=0. 

if( i-j. eq. 0)  go  to  200 

n=i-j 

val=pi*n/2. 
sval=sin(val) 
sumb=sr ( j )*sval*sval/val 
200  sum=sum+sumb 

100  s( i)=cmplx(sr( i) ,sum) 
return 
end 
c 

c  [ rossano. thesis] corr.  include 

c 

c  This  subroutine  calculates  the  correlation 


subroutine  corr  (dimt ,dimf ,s , j ,dt ,c,dp) 

integer  dimt.dimf ,dp,dp2,i, j 
complex  s(dimt) ,c(dimf) ,dum 
real  coef,dt 

dp2=dp*2 
coef=2. *dt 
do  100  i=l,dp+l 
if ( j • ge. i)  then 
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dum=s( j-i+1) 
else 

dum=cmplx(0. ,0. ) 
end  if 

c( i)=coef*(s( j+i-l)*conjg(dum)) 
if(  i.  ne.  j.  and.  i.  ne.  dp+1)  then 
c(dp2-i+2)=conjg(c( i)) 
endif 

100  continue 
return 
end 
c 

c  [ rossano. thesis] datain.  include 

c 

c  This  subroutine  reads  the  time  signal  data  into  an  array 

c 

c  ASSUME  input  file  in  format:  time  amplitude  (2x,el6. 8,5x,el6. 8) 
c  ASSUME  end  of  file  indication  is  a  last  entry  of  9999.  ,9999. 
c  Correct  format  may  be  obtained  using  [ rossano. data] convert. for 
c 


subroutine  datain  (dimt, inname, avgdelt, dp, ain,dt) 
integer  j , i ,dp, ss ,dpend,n,dimt 

real  tin(2050) ,ain(dimt) ,amp,tim,mtime,mfreq,avgdelt,df ,dt 
character*25  inname 

open  (unit=4, file=inname,status='old' ) 
rewind  4 

mt ime=dp*avgde 1 t 
df=l. / ( 4.  *mtime) 
mfreq=2.  *dp*df 

print*,'  Max  time  is' 

write(5.906)  mtime 

print*,  Delta  f  is' 

write(5j906)  df 

print*.  Max  frequency  is  now' 

write(5,906)  mfreq 

print* 

print*,'  This  can  be  changed  by  multiplying  the  delta  t  step  size' 
print*,'  if  step  size  =  1  ;  Max  freq  remains' 
write(5,906)  mfreq 

print*,'  if  step  size  =  2  ;  Max  freq  becomes  df=  mtime=' 

mtime=dp*2.  *avgdelt 

df=l. /( 4. *mtime) 

mfreq=2.  *dp*df 

write(5,907)  mfreq, df, mtime 

print*,'  If  step  size  =  5  ;  Max  freq  becomes  df=  mtime=' 

mtime=dp*5*avgdelt 

df=l. /( 4.  *mtime) 

mfreq=2. *dp*df 

write(5,907)  mfreq, df, mtime 


111 


df= 


mtime 


print*,'  If  step  size  =  10  ;  Max  freq  becomes 
mtime=dp*10*avgdelt 
df=l. /(4. *mtime) 
mfreq=2.  *dp*df 
vrite(5,907)  mf req ,df ,mtime 

print*,'  If  step  size  =  16  ;  Max  freq  becomes  df= 

mtime=dp*16*avgdelt 

df=l.  /(4. *mtime) 

mf req=2.  *dp*df 

write(5,907)  mf req ,df ,mtime 

print* 

print*,'  Input  desired  step  size' 

read(*,*)  ss 

print* 

do  111  j=l,dimt 
ain( j )=0. 0 
tin( j )=0.  0 
111  continue 

i=0 

n=ss - 1 
dpend=dp*ss 
do  200  j=l,dpend 
n=n+l 

read(4,902)  tim,amp 

if  (tim. eq. 9999. )  then 

if  ( amp.  eq.  9999.  )  go  to  210 
endif 

if  (n. eq. ss)  then 
i=i+l 

ain( i)=amp 
tin( i)=tim 
n=0 
endif 

200  continue 
go  to  212 

c  This  pads  zeros  if  the  signal  is  too  short 

210  j=i+l 

print*,'  Padding  with  zeros  since  not  enough  data  points' 
print* 

do  211  j=j ,dp 
ain( j )=0.  0 
tin( j )=j*ss*avgdelt 

211  continue 

212  close  (unit=4) 
mtime=tin( dp) 
dt=mtime/(dp-l) 

call  dtcalc  (dimt, dp, inname, tin, avgdelt) 


mt ime= 
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print*,'  Average  dt=' 
write(5.906)  avgdelt 
print*,  Final  dt=' 
write(5,906)  dt 
print* 

if  (avgdelt. ne. dt)  then 

print*,'  Average  dt  does  not  equal  Final  dt' 
print*,'  This  is  probably  because  your  time  record' 
print*,'  did  not  start  at  zero.' 

print*,'  This  can  be  fixed  by  letting  dt=average  dt' 
print*,'  Is  this  ok?' 

print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq. eq. 2)  stop 
dt=avgdelt 
endif 
return 

c  Format  statements 
902  format(2x,el6. 8,5x,el6. 8) 

906  format( f 16. 8) 

907  format( 20x, f 16. 8,fl6. 8,fl6. 8) 
end 

c 

c  [ rossano. thesis] dstaout. include 

c 

c  This  subroutine  prints  the  WDF  array  to  a  file 

c 

subroutine  dataout  (dp) 

integer  j,i,dp,dp2 
real  wdf( 1100,550) 
common  /wdfc/  wdf 

open  (unit=7 ,file='wdf.  unf ' ,status='new' , 

&form=' unformatted' ) 

print*,'  Writing  to  wdf. unf’ 

dp2=dp*2 
do  100  j=l,dp 
do  200  i=l,dp2 

write(7)  wdf(i,j) 

..00  continue 

100  continue 

close  (unit=7) 
return 
end 
c 

c  [ rossano. thesis] dataout2. include 

c 

c  This  subroutine  prints  the  RWDF  and  PWDF  arrays  to  files 

c 


subroutine  dataout2  (rwdf ,n,dp,nn,mm) 
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integer  j ,i,n,dp,dp2,nn,mm 
real  rwdf  (256,128'' 


if  (n.  eq.  2)  then 

open  (unit=7,file"'pwdf.  out1 .status='new’ ) 
print*,1  Writing  to  pwdf .  out 
endif 

if  (n. eq. 3)  then 

open  (unit=7 , f ile=’ rwdf.  out1 ,status=’new' ) 
print*,'  Writing  to  rwdf. out' 
endif 
dp2=dp*2 

do  100  j=l, dp/mm 
write(7,908)  j 
write( 7,909) 
do  200  i=l,dp2/nn 

write(7,906)  rwdf(i,j) 


200 

continue 

W’rite(  7,909) 
write(7,909) 

100 

continue 
close  (unit=7) 
return 

906 

format  (2x,el6.8) 

908 

format  (2x,i6) 

909 

format  (2x) 
end 

c 

c  [ rossano. thesis] dtcalc. include 

c 

c  This  subroutine  calculates  delta  t  from  an  input  file  (j=l)  or 

c  an  array  (j=dp) 

c 


subroutine  dtcalc  (dimt, j , inname, tin, avgdelt) 
integer  i,j,n,dimt 

real  tin(dimt) ,ain( 1025) ,sum,delt( 1025) , avgdelt 
character*25  inname 


n=0 

sum=0. 

if  ( j. ne. 1)  then 
n=j 

go  to  300 
endif 

open  (unit=4,file=inname,status='old' ) 
rewind  4 

do  100  i=l ,1025 

read(4,904)  tin( i) ,ain( i) 
if  (tin( i). eq. 9999. )  then 

if  ( ain( i). eq. 9999.  )  go  to  200 
endif 
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n=n+l 

100  continue 
200  close  (unit=4) 

300  do  400  i=l ,n-l 

delt( i)=tin( i+l)-tin( i) 
sum=sum+delt( i) 

400  continue 

avgdelt=sum/(n-l) 

return 

904  format  (2x,el6. 8,5x,el6. 8) 
end 
c 

c  [ rossano. thesis] f ft. include 

c 

c  This  subroutine  calculates  the  Fast  Fourier  Transform 

c  (FFT  if  inv=0,  Inverse  FFT  if  inv=l) 

c 

subroutine  fft  (dimf ,c,dp2,inv) 

integer  dimf , inv,dp,dp2 , i , j ,n,val,ii,k 
complex  c(dimf ) ,dum,dum2,dum3 
real  pi, const, coef 

pi=4. *atan( 1.  ) 
const=dp2 

val=alog(const)/alog(2.  )+.  1 
dp=dp2/2 

j=l 

do  40  i=l,dp2-l 

if  (i.  ge.  j)  go  to  10 
dum3=c( j ) 
c( j)=c(i) 
c( i)=dum3 
10  k=dp 

20  if  (k.  ge.  j)  go  to  30 

j=j-k 
k=k/2 
go  to  20 
30  j=j+k 

40  continue 

do  70  n=l,val 
coef=2**n 
coef=coef/2 
dum2=cmplx( 1.  ,0.  ) 

dum=cmplx( cos (pi/ coef ) , -sin(pi/coef )) 
if  (inv. ne. 0)  dum=conjg(dum) 
do  60  j=l,coef 

do  50  i=j ,dp2,2*coef 
ii=i+coef 
dum3=c( ii)*dum2 
c( ii)=c( i) -dum3 
c( i)=c( i)+dum3 
50  continue 
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dum2=dum2,vdum 

continue 

continue 

if  (inv. eq. 0)  return 
do  80  i=l,dp2 

c(i)=c(i)/cmplx( const, 0. ) 
continue 
return 
end 

[  rossano.  thesis] hamming,  include 

This  subroutine  applies  a  hamming  window  to  the  signal  ain 


subroutine  hamming  (dimt ,ain,dp,dt) 
integer  j, dp, dimt 

real  ain(dimt) ,dt ,mtime , pi , const , delta, n 

pi=4.  0*atan(  1.0) 
mtime=dp*dt 
delta=0. l*mtime 
const=pi/delta 

do  100  j=l,dp 
n=( j-l)*dt 

if  (n.  le. delta)  ain( j )=ain( j )*( .  5*( 1.  -cos( const*n) ) ) 
if  (n. ge. mtime-delta)  then 

ain( j )=ain( j )*( . 5*( 1. -cos( const*(mtime-n) ) ) ) 
endif 

100  continue 
return 
end 
c 

c  [ rossano. thesis] hanning. include 

c 

c  This  subroutine  applies  a  hanning  window  to  signal  ain 

c 

subroutine  hanning  (dimt, ain, dp, dt) 
integer  j, dp, dimt 

real  ain(dimt) ,dt,mtime,pi2,windw, const 

pi2=2.  *4.  0*atan(  1.  0) 
ratime=dp*dt 
const=(  pi2,vdt )  /mt  ime 
do  100  j=0,dp-l 

windw=0.  5*(  1  -cos(  const,vj ) ) 
ain(  j+l)=ain(  j+l),vwindw 
100  continue 
return 
end 
c 

c  [ rossano. thesis] maxamp. include 
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c 

c 

c 


This  subroutine  finds  the  maximum  amplitude  of  array  y 


subroutine  maxamp  (y,dp,zmax) 

integer  i,dp 

real  y(dp),zmax,zhold 

zmax=0.  0 

do  200  i=l,dp 
zhold=y( i) 

if  (zhold. gt.  zmax)  then 
zmax=zhold 
endif 

200  continue 

return 
end 
c 

c  [ rossano. thesis] mean,  include 

c 

c  This  subroutine  calculates  the  mean  value  in  the  signal  ain 


subroutine  mean  (dimt,ain,dp,meanv) 

integer  dp,i,dimt 
real  ain(dimt) ,sum,meanv 

sum=0. 0 
do  100  i=l,dp 
s'im=sum+ain(  i) 

100  continue 

meanv=sum/dp 

return 

end 

c 

c  [ rossano. thesis] meanr. include 

c 

c  This  subroutine  removes  the  mean  value  from  the  signal  ain 


subroutine  meanr  (dimt,ain,dp,meanv) 

integer  dp,i,dimt 
real  ain( dimt) .raeanv 

do  100  i=l,dp 

ain( i)=ain( i)-meanv 
continue 
return 
end 

[ rossano. thesis] minamp. include 
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u  o 


c 

c 


This  subroutine  finds  the  minimum  amplitude  of  array  y 


subroutine  minamp  (y, dp, zmin) 

integer  i,dp 

real  y(dp) ,zmin,zhold 

zmin=0. 0 

do  200  i=l,dp 
zhold=y( i) 

if  (zhold.  It. zmin)  then 
zmin=zhold 
end  if 

200  continue 
return 
end 

[ rossano. thesis] plot2d.  include 

c  This  subroutine  uses  disspla  11.0  to  plot  a  2  Dimension  plot 

c 

c  The  data  input  array  is  y,  the  x  axis  assumes  a  constant  dt  spacing 


subroutine  plot2d  (y,dt, dp, title, outname) 

real  dt,mtime,zmax,zmin 
character*23  datetime 
character*50  title, outname 
integer*4  status 
integer  iplot_val , dp, mesh, j 
include  1  size. include' 

call  maxamp  (y,dp,zmax) 
call  minamp  (y, dp, zmin) 

mtime=dp*dt 
do  100  j=l,dp 
x( j)=( j-l)*dt 
100  continue 

print* 

print*,  Enter  limits  for  plotting  (do  not  forget  the  decimal) 

print*,'  YMIN  YMAX' 

write(5,*)  zmin,zmax 

read(5,*)  zmin,zmax 

print* 

mesh=l 

call  reset  ( ' all' ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
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print* 

if  (iplot_val  . eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
end  if 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.’ 
call  LN03I 
endif 

call  swissm 

call  hwshd 

call  chrpat  (16) 

call  height  (.325) 

call  physor  (1.0,0.625) 

call  area2d  (6. 0,8.5) 

call  alnmes  (.5,0. ) 

call  messag  (title, 100, 3. 2,9. 9) 

call  messag  (outname, 100,3.  2,9.  5) 

call  reset  ('alnmes') 

call  blsur 

call  height  (0.200) 

call  xname  ('Time  (sec)', 10) 

call  yname  ( 'Amplitude' ,9) 

call  cross 

call  xintax 

call  yintax 

call  graf  (0.  ,mtime/4.  ,mtime,zmin,(zmax-zmin)/4.  ,zmax) 

print*,'  Axes  are  complete' 

call  reset  ('hwshd') 

call  curve  (x,y,dp,0) 

call  height  (. 150) 

call  alnmes  (0.  0,1.0) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime ,23 , -0.  8 , -0.  4) 

call  reset  ('alnmes') 

print*,'  Surface  is  complete' 

510  call  end3gr  (0) 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] plot2d2.  include 

c 

c  This  subroutine  uses  disspla  11.  0  to  plot  2  2  Dimension  curves 

c  on  one  plot 

c 

c  The  data  input  arrays  are  x  and  y,  the  x  axis  assumes  a  constant  dt 
c  spacing 
c 

subroutine  plot2d2  (x,y,dt, dp, title, outname, mnq) 
real  dt ,mtime ,zmax,zmin,xmax ,xmin,ymax ,ymin 
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character*23  datetime 
character*50  title,outname ,meani 
integer*4  status 
integer  iplot_val , dp, mesh, j ,mnq 
include  1  size. include' 


mt ime=dp*dt 
mesh=l 

if  (mnq. eq.  1)  then 

meani='Mean  value  removed' 
end  if 

if  (mnq. eq.  2)  then 

meani='Mean  value  not  removed' 
end  if 

call  maxamp  (x,dp,zmax) 
call  minamp  (x,dp,zmin) 
xmax=zmax 
xmin=zmin 

call  maxamp  (y,dp,zmax) 
call  minamp  (y,dp,zmin) 
ymax=zmax 
ymin=zmin 

if  (xmax.  gt.  ymax)  zmax=xmax 
if  (xmin.  It.  ymin)  zmin=xmin 


print* 

print*,'  Enter  limits  for  plotting  (do  not  forget  the  decimal)’ 

print*,'  YMIN  YMAX ’ 

write(5,*)  zmin.zmax 

read(5,*)  zmin.zmax 

print* 

do  100  j=l,dp 
t( j)=( j-l)*dt 
100  continue 

call  reset  ( ' all ' ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 
call  hwshd 
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call  chrpat  (16) 

call  height  ( . 325) 

call  physor  (1.0,0.625) 

call  area2d  (6. 0,8.5) 

call  alnmes  (. 5,0. ) 

call  mess ag  (title, 100,3. 2,9. 9) 

call  messag  (outname, 100,3. 2,9. 5) 

call  reset  ('alnmes') 

call  blsur 

call  height  (0.200) 

call  xname  ('Time  (sec) ',10) 

call  yname  ( 'Amplitude' ,9) 

call  cross 

call  xintax 

call  yintax 

call  graf  (0. ,mtime/4. ,mtime,zmin,(zmax-zmin)/4. ,zmax) 

print*,'  Axes  are  complete' 

call  reset  ('hwshd') 

call  curve  (t,x,dp,0) 

call  dot 

call  curve  (t,y,dp,0) 

call  height  (. 150) 

call  alnmes  (0.  0,1.0) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,-0. 8,-0. 1) 

call  messag  (meani, 50,-0. 8,-0. 4) 

call  messag  ('Solid  line=REAL' , 15,3.  0, -0. 1) 

call  messag  ('Dotted  line=IMAGINARY' ,21,3.  0,-0.  4) 

call  reset  ('alnmes') 

print*,'  Surface  is  complete' 

510  call  end3gr  (0) 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 

c 

c  [  rossano.  thesis] plot3d32.  include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  a  single  3  dimension 

c  graph  of  rwdf(64,32) 


subroutine  plot3d32  (rwdf, outname, mtime,mfreq, dp, 
&tcode,title,mnq,ampl) 

real  mtime,mfreq,zmax,zmin,ampl 
character*23 , datetime, t code, dpnum, amp 1 if 
character*50 .title , outname , meani 
integer*4  status 

integer  ip lot_va 1 , dp , mesh , rdp2 , rdp , mnq 
real  rwdf(64,32) 

rdp=32 

if  (dp. eq. 128)  then 

dpnum=’l28  data  points' 
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endif 

if  (dp.  eq.  256)  then 

dpnum-  256  data  points' 
endif 

if  (dp. eq. 512)  then 

dpnum-  512  data  points' 
endif 

if  (dp. eq. 1024)  then 

dpnum-  1024  data  points' 
endif 

if  (mnq.  eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq.  eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified' 

if  ( ampl. eq. 0. 0)  amplif='Time  amplified  by  0.0' 
if  (ampl. eq. 1. 0)  amplif='Time  amplified  by  1. 0 ' 
if  ( ampl. eq. 2. 0)  amplif='Time  amplified  by  2.0' 
if  (ampl. eq. 3. 0)  amplif='Time  amplified  by  3.0' 
if  (ampl. eq. 4. 0)  amplif='Time  amplified  by  4.0' 
if  ( ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  range32  (rwdf , rdp,zmax,zmin) 
print* 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,'  ZMIN  ZMAX' 

write(5,*)  zmin,zmax 

read(5,*)  zmin.zmax 

print* 

mesh=l 

rdp2=rdp*2 

call  reset  ('all') 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  . eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  .  eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 

call  hwshd 

call  chrpat  (16) 

call  height  (. 325) 

call  physor  (.75,. 625) 

call  area2d  (7. 5,9. 75) 

call  alnmes  (.5,0. ) 

call  messag  ( title, 100, 7. 5/2. ,9. 25) 
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call  messag  (outname, 100,7. 5/2. ,8. 75) 
call  reset  ('airlines') 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz) ',15) 

call  y3name  ('Time  (msec) ',11) 

call  z3name  (' Amplitude' ,9) 

call  vuangl  (-60.  ,30.  ,30.  ) 

call  xintax 

call  y intax 

call  zintax 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf, mesh, rdp2, mesh, rdp,0) 

call  reset  ('alnmes') 

call  height  (0. 150) 

call  alnmes  (0.  0,1.0) 

call  messag  (meani,50,5. 0,0. 0) 

call  messag  (amplif ,23,5.  0, -0.  2) 

call  mess'®  .  tcode,23,5. 0, -0. 4) 

status  =  1 ib$date_time  (datetime) 

call  eie'^ag  (datetime, 23,0.  0,0.  2) 

call  messag  (dpnum, 23,0. 0,0.  0) 

call  messag  ('Reduced  to  64  x  32 ’  ,23 , 0. 0 , -0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0. 0,-0. 4) 

call  reset  ('alnmes') 

print*,'  Surface  is  complete' 

510  call  end3gr  (0) 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano.  thesis] plot3d64.  include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  a  single  3  dimension 

c  graph  of  rwdf( 128,64) 


subroutine  plot3d64  ( rwdf .outname ,mtime ,mfreq, dp, 
&tcode,title,mnq,ampl) 

real  mtime,mfreq,zmax,zmin,ampl 
char acter*23 , datetime , tcode , dpnum , ampl if 
char act er*50 .title, outname , meani 
integer*4  status 

integer  iplot_val , dp ,mesh , rdp2 , rdp , mnq 
real  rwdf( 128,64) 

rdp=64 

if  (dp. eq. 128)  then 

dpnum='l28  data  points' 
end  if 


123 


if  (dp.  ecj.  256)  then 

dpnum-256  data  points' 
end  if 

if  (dp.  ecj.  512)  then 

dpnum-512  data  points' 
endif 

if  (dp. eq.  1024)  then 

dpnum-1024  data  points' 
endif 

if  (mnq. eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq.  eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified' 

if  ( ampl. eq. 0. 0)  amplif='Time  amplified  by  0.  O' 
if  (ampl. eq. 1. 0)  amplif='Time  amplified  by  1.  O' 
if  (ampl. eq. 2. 0)  amplif='Time  amplified  by  2.0' 
if  ( ampl. eq. 3. 0)  amplif='Time  amplified  by  3.0' 
if  ( ampl.  eq.  4.  0)  amplif^'Time  amplified  by  4.0* 
if  ( ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  range64  (rwdf , rdp,zmax,zmin) 
print* 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,'  ZMIN  ZMAX' 

write(5,*)  zmin.zmax 

read(5,*)  zmin.zmax 

print* 

mesh=l 

rdp2=rdp*2 

call  reset  ('all') 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  .  eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  sw'issm 
call  hwshd 


call  chrpat  (16) 

call  height  ( .  325) 

call  physor  (.75,.  625) 

call  area2d  (7.5,9.75) 

call  alnmes  (.5,0.  ) 

call  messag  (title, 100, 7.  5/2.  ,9.  25) 

call  messag  (outname, 100 ,7. 5/2. ,8. 75) 
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call  reset  ('alnmes') 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz)' ,15) 

call  y3name  (’Time  (msec) ’,11) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-60. ,30. ,30.  ) 

call  xintax 

call  yintax 

call  zintax 

call  graf3d  (0. ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3.  ,zmax) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf, mesh, rdp2, mesh, rdp,0) 

call  reset  ('alnmes') 

call  height  (0.150) 

call  alnmes  (0.  0,1.0) 

call  messag  (meani,50,5.  0,0.  0) 

call  messag  (amplif ,23,5. 0,-0. 2) 

call  messag  (tcode,23,5. 0,-0. 4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime,23,0. 0,0. 2) 

call  messag  ( dpnum, 23,0. 0,0. 0) 

call  messag  ('Reduced  to  128  x  64' ,23 , 0. 0 , -0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0.  0,-0. 4) 

call  reset  ('alnmes') 

print*,'  Surface  is  complete' 

510  call  end3gr  (0) 
call  endpl  (0) 
print*, '  Plotting  complete' 
print* 
return 
end 
c 

c  [  rossano.  thesis] plot3dl28. include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  a  single  3  dimension 

c  graph  of  rwdf (256, 128) 


subroutine  plot3dl28  (rwdf ,outname,mtime,mfreq,dp, 
&tcode,title,mnq,ampl) 

real  mtime,mfreq,zmax,zmin,ampl 
char act er*23 , datetime, t code, dpnum, amplif 
character*50, title, outname ,meani 
integer*4  status 

integer  iplot_val , dp ,mesh, rdp2 ,rdp ,mnq 
real  rwdf ( 256 , 128) 

rdp=128 

if  (dp. eq.  128)  then 

dpnum='l28  data  points' 
endif 

if  (dp. eq. 256)  then 
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dpnmn='256  data  points' 
end  if 

if  (dp.  eq.512)  then 

dpnum-  512  data  points' 
end  if 

if  (dp.  eq.  1024)  then 

dpnum-1024  data  points' 
endif 

if  (mnq.  eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq. eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified' 
if  ( ampl.  eq.  0.  0)  amplif= 
if  ( ampl.  eq.  1.  0)  amplif= 
if  ( ampl.  eq.  2.  0)  amplif= 
if  ( ampl.  eq.  3.  0)  amplif= 
if  ( ampl. eq.  4.  0)  amplify 
if  ( ampl. eq. 5.  0)  amplif= 


'Time  amplified  by 
'Time  amplified  by 
'Time  amplified  by 
'Time  amplified  by 
'Time  amplified  by 
'Time  amplified  by 


call  rangel28  (rwdf ,rdp,zmax,zmin) 
print* 

print*,'  Enter  limits  for  plotting  (do  not 

print*,'  ZMIN  ZMAX* 

write(5,*)  zmin,zmax 

read(5,*)  zmin,zmax 

print* 

mesh=l 


0.  O' 
1.  0’ 
2.  0’ 

3.  O' 

4.  0' 

5.  O' 


forget  decimal)' 


rdp2=rdp*2 

call  reset  ( ' all' ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy? ' 

vrite(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 

call  hwshd 

call  chrpat  (16) 

call  height  (.325) 

call  physor  (.75,. 625) 

call  area2d  (7. 5,9. 75) 

call  alnmes  (.5,0. ) 

call  messag  (title, 100, 7. 5/2.  ,9. 25) 

call  messag  (outname, 100,7. 5/2. ,8. 75) 

call  reset  ('alnmes') 
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call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz)' ,15) 

call  y3name  ('Time  (msec)', 11) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-60. ,30. ,30. ) 

cal]  xintsx 

call  yintax 

call  zintax 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 

call  reset  ('alnmes') 

call  height  (0. 150) 

call  alnmes  (0.  0,1.0) 

call  messag  (meani,50,5. 0,0.  0) 

call  messag  (amplif ,23,5. 0, -0. 2) 

call  messag  (tcode,23,5. 0,-0. 4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0. 0,0. 2) 

call  messag  ( dpnum , 23 , 0. 0 , 0. 0) 

call  messag  ('Reduced  to  256  x  128 ' , 23 , 0. 0 , -0. 2) 
call  messag  ('Smoothed  10  x  10* ,23,0. 0,-0. 4) 
call  reset  ('alnmes') 
print*,'  Surface  is  complete' 

510  call  end3gr  (0) 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] plot3dsplit32.  include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  a  series  of 

c  3  dimensional  graphs  of  rwdf(64,32) 


subroutine  plot3dsplit32  (rwdf ,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real  mtime,mfreq,zmax,zmin,rwdf( 64 , 32) , ampl 
charact er*2 3 , datetime ,t code , dpnum, amplif 
chart  cter*50 , title, outname ,meani 
integer*4  status 

integer  iplot_val ,dp ,mesh , rdp2 , rdp ,mnq 
rdp=32 

if  (dp. eq. 128)  then 

dpnum='l28  data  points' 
endif 

if  (dp. eq. 256)  then 

dpnum='256  data  points' 
endif 
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if  (dp.  eq.  512)  then 

dpnum-  512  data  points' 
end  if 

if  (dp. eq. 1024)  then 

dpnum-  1024  data  points' 
end  if 

if  (mnq.  eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq. eq. 2)  then 

meani= 'Mean  value  not  removed' 
endif 

amplif='Time  amplified' 
if  ( ampl.  eq.  0.  0)  amplif= 
if  (ampl.  eq.  1.  0)  amplif= 
if  (ampl. eq.  2.  0)  amplif= 
if  ( ampl.  eq.  3.  0)  amplif= 
if  (ampl.  eq.  4.  0)  amplif= 
if  (ampl. eq.  5.  0)  amplif= 


Time  amplified  by  0.0' 
Time  amplified  by  1.0' 
Time  amplified  by  2.0' 
Time  amplified  by  3.0' 
Time  amplified  by  4.0' 


Time  amplified  by  5.0' 
call  range32  (rwdf ,rdp,zmax,zmin) 
print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 
print*,'  ZMIN  ZMAXr 

write(5,*)  zmin,zmax 
read(5,*)  zmin,zmax 
print* 
mesh=l 
rdp2=rdp*2 
call  reset  ('all') 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy? ' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)’ 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  .  eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 
call  hwshd 
call  chrpat  (16) 


c  LABELS 

call  height  (0.200) 

call  physor  (0.5,0.625) 

call  area2d  (7. 5,9. 75) 

call  alnmes  (0. 5,0.0) 

call  raessag  (title, 100, 7. 5/2. ,9. 75) 

call  messag  (outname, 100,7.  5/2.  ,9.  5) 

call  reset  ('alnmes') 

call  height  (0. 150) 
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call  alnmes  (0.  0,1.  0) 

call  messag  (meani,50,5.  0,0.  0) 

call  messag  (amplif ,23,5. 0,-0. 2) 

call  messag  (tcode,23,5.  0,-0.  4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0.  0,0.  2) 

call  messag  ( dpnum, 23,0. 0,0. 0) 

call  messag  ('Reduced  to  64  x  32' ,23,0. 0,-0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0. 0,-0. 4) 

call  reset  ('alnmes') 

call  endgr  (1) 

print*,1  Labels  are  complete' 

c  FIRST  SUBPLOT 

call  height  (. 150) 
call  physor  (.5,5.5) 
call  area2d  (3.5,4.875) 
cal1  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz) ',15) 

call  y3name  ('Time  (msec) ',11) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-50.  ,0.  ,30.  ) 

call  xintax 

call  yintax 

call  zintax 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
ot(zmax-zmin)/3.  ,zmax) 
print*,'  First  subplot  axes  are  complete' 
call  reset  ('hwshd  ; 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 
call  endgr  (2) 

print*,'  First  subplot  surface  is  complete' 

c  SECOND  SUBPLOT 

call  height  (.150) 
call  phvsor  (4.25,5.5) 
call  area2d  (3. 5,4.5) 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ( '  '  ,  1) 

call  y3name  ('Time  (msec) ',11) 

call  z3name  (' Amplitude' , 9) 

call  vuangl  (0. ,0.  ,30.  ) 

call  yintax 

call  zintax 

call  xnonum 

call  graf3d  (0.  ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,'  Second  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 
call  endgr  (3) 

print*,'  Second  subplot  surface  is  complete' 
c  THIRD  SUBPLOT 
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call  reset  ('xnonum') 
call  height  (. 150) 
call  physor  (.5,. 625) 
call  area2d  (3.5,4.875) 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz)', 15) 

call  y3name  ( '  ' , 1) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-90. ,0.  ,30.  ) 

call  xintax 

call  zintax 

call  ynonum 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,'  Third  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf, mesh, rdp2, mesh, rdp,0) 
call  endgr  (4) 

print*,'  Third  subplot  surface  is  complete' 

C  FOURTH  SUBPLOT 

call  reset  ('ynonum') 
call  height  (. 125) 
call  physor  (4.75,1.5) 
call  area2d  (2. 75,3. 0) 
call  blsur 

call  xname  ('Frequency  (Hz)', 15) 
call  yname  ('Time  (msec) ',11) 
call  xintax 
call  yintax 

call  graf  (0.  ,mfreq/4.  ,mfreq,0.  >mtime/2.  ,mtime) 
print*,'  Fourth  subplot  axes  are  complete' 
call  messag  ('Contour  at  zmax/3' ,17,0. 5,3. 3) 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  (zmax/3.  , zmax/3.  ) 
call  conmak  (rwdf ,rdp2,rdp, 1) 
call  conend 

call  conlin  ( 1 ,' solid' , 'nolabels ', 1, 10) 
call  contur  ( 1, 'nolabels ', 'draw' ) 
call  endgr  (5) 

print*,'  Fourth  subplot  surface  is  complete' 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] plot3dsplit64. include 

c 

c  This  subroutine  uses  disspla  11.  0  to  plot  a  series  of 

c  3  dimensional  graphs  of  rwdf (128,64) 


subroutine  plot3dsplit64  ( rwdf , outname ,mtime ,mfreq ,dp, tcode , 
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&t it 1 e , mnq , amp 1 ) 

real  mtime,mfreq,zmax,zmin,rwdf( 128,64) ,ampl 
character*23 , datet ime , t code , dpnum , ampl if 
character*50 , tit le ,outname ,meani 
integer*4  status 

int  ege r  ip lot_val , dp , mesh , rdp2 , r dp , mnq 
rdp=64 

if  (dp.  ecj.  128)  then 

dpnum-  128  data  points' 
endif 

if  (dp. eq. 256)  then 

dpnum-  256  data  points’ 
endif 

if  (dp. eq.  512)  then 

dpnum-  512  data  points' 
endif 

if  (dp. eq. 1024)  then 

dpnum-1024  data  points' 
endif 

if  (mnq. eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq. eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified' 

if  ( ampl. eq. 0. 0)  amplif='Time  amplified  by  0.  O' 
if  ( ampl. eq. 1. 0)  amplif='Time  amplified  by  1.0' 
if  (ampl. eq.  2.  0)  amplif='Time  amplified  by  2.0' 
if  (  ampl.  eq.  3.  0)  amplif='Time  amplified  by  3.0' 
if  ( ampl. eq. 4. 0)  amplif='Time  amplified  by  4.0' 
if  ( ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  range64  (rwdf ,rdp,zmax,zmin) 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,'  ZMIN  ZMAX 

write(5,*)  zmin.zmax 

read(5,*)  zmin.zmax 

print* 

mesh=l 

rdp2=rdp*2 

call  reset  ( ' all'  ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 

print*,'  This  will  take  several  minutes.' 
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call  LN03I 
end  if 

call  swissm 
call  hwshd 
call  chrpat  (16) 

c  LABELS 

call  height  (0.200) 

call  physor  (0.5,0.625) 

call  area2d  (7.5,9.75) 

call  alnmes  (0. 5,0.0) 

call  messag  (title, 100,7. 5/2. ,9. 75) 

call  messag  (outname,100,7.  5/2.  ,9.  5) 

call  reset  ('alnmes') 

call  height  (0.150) 

call  alnmes  (0. 0,1.0) 

call  messag  (meani,50,5.  0,0.  0) 

call  messag  (amplif ,23,5. 0, -0. 2) 

call  messag  (tcode,23,5.  0,-0.  4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0.  0,0. 2) 

call  messag  ( dpnum, 23,0.  0,0.  0) 

call  messag  ('Reduced  to  128  x  64' ,23,0. 0, -0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0. 0,-0. 4) 

call  reset  ('alnmes') 

call  endgr  (1) 

print*,'  Labels  are  complete' 


c  FIRST  SUBPLOT 

call  height  (. 150) 
call  physor  (.5,5.5) 
call  area2d  (3.5,4.875) 
call  blsur 

call  volm3d  (8. ,8. ,9. ) 

call  x3name  ('Frequency  (Hz) ’,15) 

call  y3name  ('Time  (msec) ’,11) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-50.  ,0.  ,30.  ) 

call  xintax 

call  yintax 

call  zintax 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3.  ,zmax) 
print*,'  First  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf, mesh, rdp2, mesh, rdp,0) 
call  endgr  (2) 

print*,'  First  subplot  surface  is  complete' 


c  SECOND  SUBPLOT 
call  height 
call  physor 
call  area2d 
call  blsur 
call  volm3d 
call  x3name 


(• 150) 

(4.25.5.5) 

(3.5. 4. 5) 

(8.  ,8.  ,9.  ) 

(’  M) 
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call  y3name  ('Time  (msec)', 11) 
call  z3name  ( 'Amplitude' ,9) 
call  vuangl  (0.  ,0.  ,30.  ) 
call  yintax 
call  zintax 
call  xnonum 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,otime,zmin, 
&(zmax-zmin)/3.  ,zmax) 
print*,'  Second  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf, mesh, rdp2, mesh, rdp,0) 
call  endgr  (3) 

print*,'  Second  subplot  surface  is  complete' 


c  THIRD  SUBPLOT 

call  reset  ('xnonum') 
call  height  (. 150) 
call  physor  (.5,. 625) 
call  area2d  (3.5,4.875) 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz) ’,15) 

call  y3name  ( '  '  ,  1) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-90.  ,0.  ,30.  ) 

call  xintax 

call  zintax 

call  ynonum 

call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3.  ,zmax) 
print*,'  Third  subplot  axes  are  complete' 
call  reset  (’hwshd  j 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 
call  endgr  (4) 

print*,'  Third  subplot  surface  is  complete' 


C  FOURTH  SUBPLOT 

call  reset  ('ynonum') 
call  height  ( . 125) 
call  physor  (4.75,1.5) 
call  area2d  (2.75,3.0) 
call  blsur 

call  xname  ('Frequency  (Hz) ',15) 
call  yname  ('Time  (msec) ',11) 
call  xintax 
call  yintax 

call  graf  (0.  ,mfreq/4.  ,mfreq,0. ,mtime/2.  ,mtime) 
print*,'  Fourth  subplot  axes  are  complete' 
call  messag  ('Contour  at  zmax/3' ,17,0. 5,3. 3) 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  (zmax/3.  , zmax/3.  ) 
call  conraak  (rwdf ,rdp2,rdp,l) 
call  conend 

call  conlin  ( 1 ,' solid ', 'nolabels 1 , 10) 
call  contur  ( 1, 'nolabels' , 'draw' ) 
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call  endgr  (5) 

print*,'  Fourth  subplot  surface  is  complete' 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] plot3dsplitl28. include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  a  series  of 

c  3  dimensional  graphs  of  rwdf(256, 128) 


subroutine  plot3dsplitl28  (rwdf ,outname,mtime,mfreq,dp,tcode, 
&t it 1 e , mnq , amp 1 ) 

real  mtime,mfreq,zmax,zmin, rwdf (256, 128) ,ampl 
char acter*23 , datetime, tcode , dpnum, ampl if 
character*50 .title, outname , meani 
integer*4  status 

integer  iplot_val , dp , mesh , rdp2 , rdp , mnq 
rdp=128 

if  (dp.  eq.  128)  then 

dpnum-  128  data  points' 
end  if 

if  (dp.  ecj.  256)  then 

dpnum-  256  data  points' 
end  if 

if  (dp. eq.  512)  then 

dpnum='512  data  points' 
end  if 

if  (dp. eq. 1024)  then 

dpnum-  1024  data  points' 
end  if 

if  (mnq.  eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq.  eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified' 

if  (ampl. eq. 0. 0)  amplif='Time  amplified  by  0.  O' 
if  (ampl. eq. 1. 0)  amplif='Time  amplified  by  1.0' 
if  (ampl. eq. 2. 0)  amplif='Time  amplified  by  2.0' 
if  (ampl. eq. 3. 0)  amplif='Time  amplified  by  3.0' 
if  ( ampl. eq. 4. 0)  amplif='Time  amplified  by  4.0' 
if  ( ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  rangel28  (rwdf , rdp, zmax,zmin) 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,'  ZMIN  ZMAXr 

write(5,*)  zmin,zmax 

read(5,*)  zmin,zmax 

print* 

mesh=l 
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rdp2=rdp*2 

call  reset  ( 'all' ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  .  eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 
call  hwshd 
call  chrpat  (16) 


c  LABELS 

call 
call 
call 
call 
call 
call 
call 
call 
call 
call 
call 
call 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0. 0,0. 2) 

call  messag  (dpnum,23,0. 0,0. 0) 

call  messag  ('Reduced  to  256  x  128' ,23,0. 0, -0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0.  0,-0. 4) 

call  reset  ('alnmes') 

call  endgr  (1) 

print*,'  Labels  are  complete' 


height  (0.200) 

physor  (0. 5,0. 625) 

area2d  (7. 5,9. 75) 

alnmes  (0. 5,0.0) 

messag  ( tit le, 100 , 7. 5/2. ,9. 75) 

messag  (outname, 100,7.  5/2. ,9. 5) 

reset  ( ' alnmes' ) 

height  (0. 150) 

alnmes  (0. 0,1. 0) 

messag  (meani ,50 ,5. 0 ,0.  0) 

messag  ( amplif ,23 ,5.  0 , -0.  2) 

messag  ( tcode , 23 ,5. 0 , -0. 4) 


c  FIRST  SUBPLOT 

call  height 
call  physor 
call  area2d 
call  blsur 
call  volm3d 
call  x3name 
call  y3name 
call  z3name 
call  vuangl 
call  xintax 
call  yintax 
call  zintax 


(• 150) 

(.5,5.5) 

(3.5,4.875) 

(8.  ,8.  ,9.) 

('Frequency  (Hz)', 13) 
( 'Time  (msec) ' , 11) 

( ' Amplitude' ,9) 

(-50.  ,0.  ,30.  ) 
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call  graf3d  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,'  First  subplot  axes  are  complete' 
call  reset  ('hwshd  j 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 
call  endgr  (2) 

print*,'  First  subplot  surface  is  complete' 

c  SECOND  SUBPLOT 

call  height  (. 150) 
call  physor  (4.25,5.5) 
call  area2d  (3.5,4. 5) 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ( '  ,1) 

call  y3name  ('Time  (msec) ’,11) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (0.  ,0.  ,30.  ) 

call  yintax 

call  zintax 

call  xnonum 

call  graf3d  (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&( zmax-zmin)/3. ,zmax) 
print*,'  Second  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  ( rwdf , mesh, rdp2 , mesh, rdp,0) 
call  endgr  (3) 

print*,'  Second  subplot  surface  is  complete' 

c  THIRD  SUBPLOT 

call  reset  ('xnonum') 
call  height  (. 150) 
call  physor  (.5,. 625) 
call  area2d  (3.5,4.875) 
call  blsur 

call  volm3d  (8.  ,8.  ,9.  ) 

call  x3name  ('Frequency  (Hz) ',15) 

call  y3name  ( '  ' , 1) 

call  z3name  ( 'Amplitude' ,9) 

call  vuangl  (-90. ,0. ,30. ) 

call  xintax 

call  zintax 

call  ynonum 

call  graf3d  (0. ,mfreq/4.  ,mfreq,0. ,mtime/2.  ,mtime,zmin, 
&(zmax-zmin)/3.  ,zmax) 
print*,'  Third  subplot  axes  are  complete' 
call  reset  ('hwshd') 

call  surmat  (rwdf , mesh, rdp2, mesh, rdp,0) 
call  endgr  (4) 

print*,'  Third  subplot  surface  is  complete' 

C  FOURTH  SUBPLOT 

call  reset  ('ynonum') 
call  height  ( . 125) 
call  physor  (4.75,1.5) 
call  area2d  (2. 75,3. 0) 
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call  blsur 

call  xname  ('Frequency  (Hz) ’,15) 
call  yname  ('Time  (msec) ',11) 
call  xintax 
call  yintax 

call  graf  (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime) 
print*,'  Fourth  subplot  axes  are  complete' 
call  messag  ('Contour  at  zmax/3' , 17 ,0. 5 ,3. 3) 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  (zmax/3.  , zmax/3.  ) 
call  conmak  (rwdf ,rdp2,rdp, 1) 
call  conend 

call  conlin  ( 1 ,' solid' , 'nolabels ', 1 , 10) 
call  contur  ( 1 , 'nolabels ',' draw' ) 
call  endgr  (5) 

print*,'  Fourth  subplot  surface  is  complete' 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano.  thesis] plotcon32.  include 

c 

c  This  subroutine  uses  Sisspla  11.0  to  plot  contour 

c  graphs  of  rwdf(64,32) 


subroutine  plotcon32  (rwdf ,outname,mtime,mfreq,dp,tcode, 
&title, mnq, ampl) 

real  mtime,mfreq,zmax,zmin,rwdf(64,32) ,ampl 
char act er*23 ,datet ime,t code ,dpnum, ampl if 
character*50 , title ,outname,meani 
integer*4  status 

integer  iplot_val ,dp,mesh,rdp2,rdp,mnq 
rdp=32 

if  (dp. eq. 128)  then 

dpnum='l28  data  points' 
endif 

if  (dp. eq. 256)  then 

dpnum=’256  data  points' 
endif 

if  (dp. eq.  512)  then 

dpnum='512  data  points' 
endif 

if  (dp. eq.  1024)  then 

dpnum='l024  data  points' 
endif 

if  (mnq. eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq. eq. 2)  then 

meani='.M€an  value  not  removed' 
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endif 

amplif='Time  amplified' 

if  ( ampl. eq. 0. 0)  amplif=' Time  amplified  by  0. O' 
if  ( ampl. eq.  1.  0)  amplif='Time  amplified  by  1.0' 
if  ( ampl.  eq.  2.  0)  amplif='Time  amplified  by  2.0' 
if  ( ampl. eq. 3. 0)  amplif='Time  amplified  by  3.0' 
if  ( ampl. eq. 4. 0)  amplif=' Time  amplified  by  4.0' 
if  (ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  range32  ( rwdf ,rdp,zmax,zmin) 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,'  ZMIN  ZMAXr 

write(5,*)  zmin,zmax 

read(5,*)  zmin,zmax 

print* 

mesh=l 

rdp2=rdp*2 

call  reset  ('all') 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  . eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 
call  hwshd 
call  chrpat  (16) 

c  LABELS 

call  height  (0.200) 

call  physor  (0.5,0.625) 

call  area2d  (7.5,9.75) 

call  alnmes  (0. 5,0. 0) 

call  messag  (title, 100, 7. 5/2. ,9. 75) 

call  messag  (outname, 100,7. 5/2. ,9. 5) 

call  messag  ('Contours  from  zmax/6  to  zmax' , 100, 7. 5/2. ,9. 25) 

call  reset  ('alnmes') 

call  height  (0. 150) 

call  alnmes  (0. 0,1.0) 

call  messag  (meani,50,5. 0,0. 0) 

call  messag  (amplif ,23,5. 0,-0. 2) 

call  messag  (tcode,23,5. 0,-0. 4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0. 0,0. 2) 

call  messag  ( dpnum, 23,0. 0,0. 0) 

call  messag  ('Reduced  to  64  x  32' ,23,0. 0,-0. 2) 

call  messag  ('Smoothed  10  x  10 ' ,23,0. 0,-0. 4) 

call  reset  ('alnmes') 
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call  endgr  (1) 

print*,'  Labels  are  complete' 


c  PLOT 

call  height  (. 125) 
call  physor  (0.75,1.5) 
call  area2d  (7.0,8.  G) 
call  blsur 

call  xname  ('Frequency  (Hz) ',15) 

call  yname  ('Time  (msec) ',11) 

call  xintax 

call  yintax 

call  xticks  (10) 

call  yticks  (10) 

call  graf  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  (zmax/6.  ,zmax) 

call  conmak  (rwdf ,rdp2,rdp, zmax/6.  ) 

call  conend 

call  conlin  ( 1 ,' solid' ,' nolabels ', 1 , 10) 
call  conlin  (2 ,' dot ',' nolabels ', 1 , 10) 
call  conlin  ( 3 ,' solid' ,' nolabels ', 1 , 10) 
call  conlin  (4, 'dot' , 'nolabels' ,1,10) 
call  conlin  (5, ' solid’ , 'nolabels' , 1 , 10) 
call  conlin  ( 6, ' dot ', 'nolabels ' ^ 1 , 10) 
call  contur  (6, 'nolabels' , 'draw' ) 
call  endgr  (2) 

print*,'  Contours  are  complete' 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] plotcon64.  include 

c 

c  This  subroutine  uses  disspla  11.0  to  plot  contour 

c  graphs  of  rwdf (128,64) 


subroutine  plotcon64  (rwdf ,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real  mtime,mfreq,zmax,zmin,rwdf( 128,64) ,ampl 
character*23, date time, tcode ,dpnum, amplif 
character*50 , title ,outname ,meani 
integer*4  status 

integer  ip lot_va 1 , dp , mesh , rdp2 , r dp , mnq 
rdp=64 

if  (dp.  eq.  128)  then 

dpnum-  128  data  points’ 
endif 

if  (dp.  eq.  256)  then 
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dpnum='256  data  points' 
end  if 

if  (dp. eq. 512)  then 

dpnum-  512  data  points' 
endif 

if  (dp.  ecp  1024)  then 

dpnum-1024  data  points' 
endif 

if  (mnq.  eq.  1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq. eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

arr.plif='Time  amplified' 
if  ( ampl. eq.  0.  0)  amplif= 

( ampl.  eq.  1.0) 

( ampl.  eq.  2.0) 

( ampl.  eq.  3.0) 

( ampl.  eq.  4.0) 

( ampl.  eq.  5.0) 


'Time  amplified  by  0.  O' 

'Time  amplified  by  1.0' 

'Time  amplified  by  2.0' 

'Time  amplified  by  3.0' 

'Time  amplified  by  4.0' 

:'Time  amplified  by  5.0' 
call  range64  (rwdf ,rdp,zmax,zmin) 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal) 
'  7MTV  ZMAX' 


if 

if 

if 

if 

if 


amplif= 
amplif= 
amplif= 
amplif= 
ampl if =' 


print*,'  ZMIN 

write(5,*)  zmin,zmax 
read(5,*)  zmin,zmax 
print* 
mesh=l 


rdp2=rdp*2 

call  reset  ('all') 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or 
&  a  hardcopy? ' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  . eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 
print*,'  This  will  take  several  minutes.' 
call  LN03I 
endif 

call  swissm 
call  hwshd 
call  chrpat  (16) 


get 


c  LABELS 

call 

call 

call 

call 

call 

call 


height 

(0.  200) 

physor 

(0.  5,0.  625) 

area2d 

(7.5,9.  75) 

alnmes 

(0.  5,0.  0) 

messag 

(title, 100, 7. 5/2.  ,9. 

75) 

messag 

(outname, 100,7.  5/2.  , 

9.5) 
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call  messag  ('Contours  from  zmax/6  to  zmax' , 100,7. 5/2. ,9. 25) 

call  reset  ('alnmes') 

call  height  (0. 150) 

call  alnmes  (0.  0,1.0) 

call  messag  (meani,50,5. 0,0. 0) 

call  messag  ( amplif , 23 ,5. 0 , -0.  2) 

call  messag  (tcode,23,5. 0,-0. 4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0.  0,0.  2) 

call  messag  (dpnum, 23,0. 0,0.  0) 

call  messag  ('Reduced  to  128  x  64' ,23,0. 0, -0. 2) 

call  messag  ('Smoothed  10  x  10' ,23,0. 0,-0. 4) 

call  reset  ('alnmes') 

call  endgr  (1) 

print*,'  Labels  are  complete' 


c  PLOT 

call  height  ( . 125) 
call  physor  (0.75,1.5) 
call  area2d  (7.0,8. 0) 
call  blsur 

call  xname  ('Frequency  (Hz) ',15) 

call  yname  ('Time  (msec)', 11) 

call  xintax 

call  yintax 

call  xticks  (10) 

call  yticks  (10) 

call  graf  (0.  ,mfreq/4.  ,mfreq,0.  ,mtime/2.  ,mtime) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  i zmax/6. ,zmax) 

call  conmak  (rwdf ,rdp2,rdp, zmax/6.  ) 

call  conend 

call  conlin  ( 1 ,' solid' ,' nolabels ', 1 , 10) 
call  conlin  (2, 'dot' , 'nolabels ', 1 , 10) 
call  conlin  ( 3 ,' solid' ,' nolabels ', 1 , 10) 
call  conlin  (4, ' dot ',' nolabels ', 1 , 10) 
call  conlin  ( 5 ,' solid' , 'nolabels ', 1 , 10) 
call  conlin  (6, ' dot' , 'nolabels' , 1 , 10) 
call  contur  (6, 'nolabels' , 'draw' ) 
call  endgr  (2) 

print*,'  Contours  are  complete' 
call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano.  thesis] plotconl28.  include 

c 

c  This  subroutine  uses  disspla  11. 0  to  plot  contour 

c  graphs  of  rwdf(256, 128) 


subroutine  plotconl28  (rwdf ,outname,mtime,mfreq,dp,tcode, 
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&t it le , mnq , amp 1 ) 

real  mtime,mfreq,zmax,zmin,rwdf(256, 128) ,ampl 
char acter*23,datetime,tcode,dpnum, ampl if 
character*50,title,outname,meani 
integer*4  status 

integer  iplot_val , dp, mesh , rdp2 , rdp ,mnq 
rdp=128 

if  (dp.  eq.  128)  then 

dpnum-  128  data  points' 
end  if 

if  (dp.  ecj.  256)  then 

dpnum=  256  data  points' 
endif 

if  (dp. eq. 512)  then 

dpnum-  512  data  points' 
endif 

if  (dp. eq. 1024)  then 

dpnum-  1024  data  points' 
endif 

if  (mnq. eq. 1)  then 

meani='Mean  value  removed' 
endif 

if  (mnq.  eq.  2)  then 

meani='Mean  value  not  removed' 
endif 

amplif='Time  amplified’ 

if  ( ampl. eq. 0. 0)  amplif='Time  amplified  by  0.  O' 
if  ( ampl. eq. 1. 0)  amplif='Time  amplified  by  1 .  0 ' 
if  ( ampl. eq. 2. 0)  amplif='Time  amplified  by  2.0' 
if  (ampl. eq. 3.  0)  amplif='Time  amplified  by  3.0' 
if  ( ampl. eq. 4. 0)  amplif='Time  amplified  by  4.0' 
if  (ampl. eq. 5. 0)  amplif='Time  amplified  by  5.0' 
call  rangel28  (rwdf ,rdp,zmax,zmin) 

print*,'  Enter  limits  for  plotting  (do  not  forget  decimal)' 

print*,’  ZMIN  ZMAX' 

write(5,*)  zmin,zmax 

read(5,*)  zmin,zmax 

print* 

mesh=l 

rdp2=rdp*2 

call  reset  ( ' all'  ) 

write(5,*)  '  Do  you  want  to  view  the  plot  on  the  screen  or  get 
&  a  hardcopy?' 

write(5,*)  '  (1  for  view,  2  for  hardcopy)' 

print* 

read(5,*)  iplot_val 
print* 

if  (iplot_val  .  eq.  1)  then 

print*,'  When  finished  viewing  hit  return  key' 
call  pgpx 
endif 

if  (iplot_val  . eq.  2)  then 

print*,'  Please  be  patient' 

print*,'  This  will  take  several  minutes.' 
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call  LN03I 
endii 

call  swissm 
call  hwshd 
call  chrpat  (16) 


c  LABELS 

call  height  (0.200) 

call  physor  (0.5,0.625) 

call  area2d  (7.5,9.75) 

call  alnmes  (0. 5,0. 0) 

call  messag  (title , 100 , 7.  5/2.  , 9. 75) 

call  messag  (outname,100,7.  5/2.  ,9.5) 

call  messag  ('Contours  from  zmax/6  to  zmax' , 100,7.  5/2. 

call  reset  ('alnmes') 

call  height  (0. 150) 

call  alnmes  (0.  0,1.0) 

call  messag  (meani,50,5.  0,0.  0) 

call  messag  (amplif ,23,5.  0,-0.  2) 

call  messag  (tcode,23,5.  0,-0.  4) 

status  =  lib$date_time  (datetime) 

call  messag  (datetime, 23,0.  0,0.  2) 

call  messag  (dpnum, 23,0.  0,0.  0) 

call  messag  ('Reduced  to  256  x  128' ,23,0. 0,-0. 2) 
call  messag  ('Smoothed  10  x  10' ,23,0. 0,-0. 4) 
call  reset  ('alnmes') 
call  endgr  (1) 

print*,'  Labels  are  complete' 


c  PLOT 

call  height  (. 125) 
call  physor  (0.75,1.5) 
ca)l  area2d  (7.0,8.  0) 
call  blsur 

call  xname  ('Frequency  (Hz) ',15) 

call  yname  ('Time  (msec)', 11) 

call  xintax 

call  yintax 

call  xticks  (10) 

call  yticks  (10) 

call  graf  (0.  ,mfreq/4.  ,mfreq,0.  ,ratime/2.  ,mtime) 
print*,'  Axes  are  complete' 
call  reset  ('hwshd') 
call  conbgn 

call  zrange  (zmax/6. ,zmax) 

call  conmak  (rwdf ,rdp2,rdp, zmax/6.  ) 

call  conend 

call  conlin  ( 1, 'solid' , 'nolabels ', 1, 10) 
call  conlin  ( 2 ,' dot ', 'nolabels ', 1 , 10) 
call  conlin  (3, ' solid' , 'nolabels ', 1 , 10) 
call  conlin  (4, ' dot' , ’nolabels' , 1 , 10) 
call  conlin  ( 5 ,' solid' , 'nolabels ', 1 , 10) 
call  conlin  (6 ,' dot ', 'nolabels ', 1 , 10) 
call  contur  ( 6, 'nolabels' ,' draw' ) 
call  endgr  (2) 

print*,'  Contours  are  complete' 


.25) 
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call  endpl  (0) 
print*,'  Plotting  complete' 
print* 
return 
end 
c 

c  [ rossano. thesis] pseudo,  include 

c 

c  This  subroutine  smooths  the  WDF  along  both  the  time  and 

c  frequency  axes 

c 


subroutine  pseudo  (dimt,dimf ,rwdf ,dp,dt,nn,mm) 

integer  dp,i, j ,dimt ,dimf ,dp2,rdp,rdp2,nn,mm,nf ,mt ,nf2,mt2, 
&ii, j j ,L,LL,k,kk 

real  rwdf(256, 128) ,hg( -25:  25,-25:  25) ,pi,dt,df , 

&wdf( 1100, 550), coef.val 
common  /wdfc/  wdf 

pi=4. *atan( 1.  ) 
df=l.  /(4.  *dp*dt) 
dp2=dp*2 
rdp=dp/mm 
rdp2=dp2/nn 

print*,'  Smoothing  10  x  10' 

nf=10 

mt=10 

nf2=nf*2 

mt2=mt*2 

val=l. / ( C  2. *pi*nf*df )*(2. *pi*mt*dt)) 
do  20  j=-mt2,mt2 
do  10  i=-nf2,nf2 

coef=-( ( j*j )/( 2. *mt*mt) )-(( i* i ) / ( 2 . *nf*nf ) ) 
hg( i , j )=val*exp( coef ) 


10 

continue 

20 

continue 

do  100  j=l,rdp 

do  100  i=l,rdp2 

100 

rwdf( i, j)=0.  0 

do  500  i=l,dp2,nn 

ii=l+( i-l)/nn 


do  450  j=l,dp,mm 
jj=l+( j-l)/mm 
do  400  L=i-nf2 , i+nf2 
LL=L 

if  (L.  It.  1)  LL=l+dp2 
if  (L.  gt.  dp2)  LIr=l-dp2 
do  350  k=j-rat2, j+mt2 
kk=k 

if  (k. It. 1)  kk=k+dp 
if  (k. gt. dp)  kk=k-dp 

rwdf(ii, j j)=rwdf(ii, j j)+wdf(LL,kk)*hg(L-i,k-j) 
350  continue 

400  continue 

450  continue 
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500  continue 

c  This  loop  removes  a  plotting  artifact 
do  600  j=l,rdp 

rwdf(rdp2-3, j)=0.  0 
rwdf ( rdp2-2 , j )=0.  0 
rwdf(rdp2-l, j)=0.  0 
rwdf(rdp2 , j )=0.  0 
600  continue 
return 
end 
c 

c  [ rossano. thesis] range32.  include 

c 

c  This  subroutine  finds  the  maximum  and  minimum  amplitudes 

c  of  array  rwdf(64,32) 


subroutine  range32  (rwdf , dp, zmax, zmin) 

real  zmax, zmin, rwdf(64, 32) 
integer  i,j,dp,dp2 

dp2=dp*2 
zmax=0.  0 
zmin=0.  0 
do  100  j=l,dp 
do  200  i=l ,dp2 

if  ( rwdf ( i,j ). gt. zmax)  zmax=rwdf(i, j) 
if  (rwdf(i,  j).  It.  zmin)  znin=rwdf ( i , j ) 

200  continue 
100  continue 
return 
end 
c 

c  [  rossano.  thesis] range64.  include 

c 

c  This  subroutine  fin^-.  the  maximum  and  minimum  amplitudes 

c  of  array  rwdf (128,6 


subroutine  range64  (rwdf , dp, zmax, zmin) 

real  zmax, zmin, rwdf(128, 64) 
integer  i,j,dp,dp2 

dp2=dp*2 
zmax=0.  0 
zmin=0. 0 
do  100  j=l,dp 
do  200  i=l,dp2 

if  (rwdf( i, j).  gt.  zmax)  zmax=rwdf ( i , j ) 
if  (rwdf( i, j ).  It.  zmin)  zmin=rwdf( i, j ) 
200  continue 
100  continue 
return 
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end 


c  [ rossano. thesis] rangel28. include 

c 

c  This  subroutine  finds  the  maximum  and  minimum  amplitudes 

c  of  array  rwdf (256 , 128) 


subroutine  rangel28  (rwdf ,dp,zmax,zmin) 

real  zmax,zmin, rwdf (256, 128) 
integer  i,j,dp,dp2 

dp2=dp*2 
zmax=0.  0 
zmin=0.  0 
do  100  j=l,dp 
do  200  i=l ,dp2 

if  (rwdf(i,  j).  gt.  zrnax)  zmax=rwdf(i, j) 
if  (rwdf(i, j).  It.  zmin)  zrain=rwdf( i, j ) 

200  continue 

100  continue 
return 
end 
c 

c  [ rossano. thesis] reduce. include 

c 

c  This  subroutine  reduces  the  data  in  wdf  to  rwdf 


subroutine  reduce  (dimt ,dimf ,dt ,dp ,nn,mm,rwdf ,rval) 

integer  dimt ,dimf ,dp,dp2 , i , j , ii, j j ,nn,mm,rval 
real  wdf ( 1100,550) , rwdf (256 , 128) 
common  /wdfc/  wdf 

dp2»dp*2 
df=l.  /(4.  *dp*dt) 


print*, ' 

What  do 

you 

want  to  reduce  to?' 

print*, ' 

1. 

64 

x  32' 

print*, ' 

2. 

128 

x  64' 

print*, ' 

3. 

256 

x  128' 

read(5,*)  rval 

if  (rval.eq.  1)  then 
nn=dp2/64 
mm=dp/32 
endif 

if  (rval.eq.  2)  then 
nn=dp2/128 
mm=dp/64 
endif 

if  (rval.eq.  3)  then 
nn=dp2/256 
mm=dp/ 128 
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endif 

ii=0 

jj=0 

do  100  j = 1 , dp , mm 
ii=0 

do  200  i=l,dp2,nn 
ii=ii+l 

rwdf(ii, j j)=wdf(i, j) 

200  continue 
100  continue 
return 
end 
c 

c  [  rossano. thesis] size,  include 

c 

c  This  code  is  included  in  wigfunl. for  and  wigfunlb.  for  and  their 

c  subroutines.  It  provides  an  easy  way  to  change  the  size  of  the 

c  time  plotting  arrays, 

c 

c  real  x( 128) ,y( 128) ,t( 128) 

c  real  x(256) ,y(256) ,t(256) 

real  x(512) ,y(512) ,t(512) 
c  real  x(1024) ,y( 1024) ,t( 1024) 

integer  sizedp 

sizedp=512 

c 

c  [ rossano. thesis] timesig.  include 

c 

c  This  subroutine  modifies  and  plots  the  time  signal 

c 


subroutine  timesig  (dimt , ain, dp, dt, outname, tcode, 
&title,s,sr,si, mnq , amp 1 ) 

integer  dp , mnq .wndwcode , anq ,  dimt , atvq 
real  ain(dimt) ,sr(dimt) ,si(dirat) ,meanv,dt ,zmax,ampl 
complex  s(dimt) 
character*23  tcode 

character*50  tit le , tit lehold , outname 
print* 

print*,'  Do  you  want  to  plot  the  raw  time  signal?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq.  eq.  1)  then 

title='Raw  Time  Signal' 
call  plot2d  (ain,dt, dp, title, outname) 
endif 
print* 

call  mean  (dimt ,ain,dp,meanv) 
print*,'  The  mean  value  is' 
write(5,906)  meanv 
print* 
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print*,’  Do  you  want  to  remove  the  mean  from  the  time  signal?' 
print*,’  (1  for  yes,  2  for  no)' 

read(*,*)  mnq 

if  (mnq. eq. 1)  call  meanr  (dimt , ain,dp,meanv) 

call  maxamp  (ain,dp,zmax) 

print* 

print*,’  The  max  amplitude  is' 

print*,zmax 

print* 

print*,'  You  want  this  to  be  approx  1  in  order  to  avoid 
&plotting  artifacts.  ' 
ampl=l.  0/zmax 
print* 

print*,'  Recommend  an  amplification  of' 

print*,ampl 

print* 

print*,'  Do  you  want  to  amplify  the  signal?' 
print*,'  (1  for  yes,  2  for  no)’ 

read(5,*)  atvq 

if  (atvq. eq. 1)  call  amplify  (dimt,ain,dp,ampl) 

if  (atvq.  ne.  1)  amp  1=0.  0 

print* 

print*,'  Do  you  want  to  window  the  time  signal?' 
print*,'  (1  for  yes,  2  for  no)r 

read(5,*)  atvq 
if  (atvq.  eq. 1)  then 

call  window  (dimt,ain,dp,dt,wndwcode) 
if  (wndwcode.  eq.  2)  tcode=' Hanning  window  time' 
if  (wndwcode.  eq.  3)  tcode=l Hamming  window’  time’ 
endif 
print* 

print*,'  Do  you  want  to  plot  the  modified  time  signal?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq. eq. 1)  then 

title='Modified  Time  Signal' 
call  plot2d  ( ain.dt , dp, title, outname) 
endif 

title=' Wigner  Distribution' 
print* 

print*,'  Do  you  want  to  make  the  time  signal  analytic?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  anq 
if  (anq. eq. 1)  then 

title='Wigner-Ville  Distribution' 
call  analytic  (dimt ,ain,s ,dp) 
do  100  j=l,dp 
sr( j)=ain( j) 
si( j )=aimag( s( j ) ) 

100  continue 
print* 

print*,'  Do  you  want  to  plot  the  analytic  time  signal?' 
print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq.  eq. 1)  then 
print* 
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if  (dp.  ge.  512)  then 

print*,'  The  STD00001.dat  file  may  be  too  large  to  print.  ' 
print*,'  For  a  hardcopy  run  [ rossano.  thesis] wigfunlb.  for ' 
print* 

print*,'  Do  you  want  to  continue  with  analytic  plotting 

&  here?' 

print*,'  (1  for  yes,  2  for  no)' 

read(5,*)  atvq 
if  (atvq. eq. 2)  go  to  300 
endif 

titlehold=title 
title=' Analytic  Time  Signal' 
call  plot2d2  (sr ,si,dt ,dp,title,outname,mnq) 
title=titlehold 
endif 
endif 

300  if  (anq. eq. 2)  then 
do  200  j=l,dp 
sr( j)=ain( j) 
si(  j)=0.  0 

s( j)=cmplx(sr( j)  ,si( j)) 


200 

continue 

endif 

return 

906 

format( f 16. 8) 

end 

c  [ rossano.  thesis] wigh.  include 

c 

c  This  subroutine  calculates  the  WDF 


subroutine  wigh  (dimt,dimf ,s,c,df ,dt,dp2,dp) 

integer  dimt ,dimf , i , j ,dp2 ,dp 
complex  s(dimt) ,c(dimf ) 
real  wdf( 1100,550) ,df,dt 
common  /wdfc/  wdf 

do  100  j=l,dp 

call  corr  (dimt,dimf ,s , j ,dt,c,dp) 
call  fft  (dimf ,c,dp2,0) 
do  200  i=l,dp2 

wdf ( i , j )=real( c( i) ) 

200  continue 

100  continue 
return 
end 
c 

c  [ rossano. thesis] window. include 

c 

c  This  subroutine  calls  the  available  windowing  functions 


subroutine  window  (dimt,ain,dp,dt,wndwcode) 
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integer  atvq,dp,wndwcode,dimt 
real  ain(dimt) ,dt 

print* 

print*,'  Which  window  would  you  like  to  apply?' 

print*,'  1.  None' 

print*,'  2.  Hanning' 

print*,'  3.  Hamming' 

read(5,*)  atvq 

if  (atvq. eq. 1)  go  to  100 

if  (atvq. eq. 2)  then 

call  hanning  (dimt,ain,dp,dt) 
wndwcode=2 
end  if 

if  (atvq. eq.  3)  then 

call  hamming  (dimt,ain,dp,dt) 
wndwcode=3 
end  if 

100  return 
end 

D.  DATA  CONVERSION  PROGRAM 

c 

c  [ rossano. thesis] convert,  for 

c 

c  This  program  converts  the  experimental  input  data  file  obtained 

c  from  HP  Vista  into  the  format  used  for  WIGFUN1.F0R 


real  tim(2) ,amp(2) 
integer  modq,i,n 
character*20  inname, outname 
character*l  c(6) 

print* 

print*,'  This  program  converts  HP  Vista  data  input  format  into' 
print*,'  the  data  format  used  in 
&  [ rossano. thesis] VIGFUN1.  FOR' 
print* 

print*,'  Enter  HP  Vista  input  filename' 

read(5,901)  inname 

print* 

print*,'  Has  the  last  line  been  modified  so  that  9.999999' 
print*,'  is  in  each  column  to  indicate  an  end  of  file?' 
print*,'  An  example  follows: ' 
print*,'  (  4. 960938e-01,  7. 668274e-01) 

&  (  4. 980469e-01 ,  1. 742336e-01) ' 

print*,'  (  5. 000000e-01,  -5. 128824e-01) 

&  (  9.999999  ,  9.999999  )' 

print*,’  (  9.999999  ,  9.999999  ) 

&  (  9.999999  ,  9.999999  )' 

print*,'  (1  for  yes,  2  for  no)' 

read( 5 ,902)modq 
if  (modq. eq.  2)  then 

print*,'  You  need  to  edit  the  HP  Vista  file’ 
print*,’  Good  bye' 
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go  to  900 
endif 
print* 

print*,’  Enter  output  filename  (i. e.  sin62.5e)' 
read(5,901)  outname 

open  (unit=4,file=inname,status='old' ) 
rewind  4 

open  (unit=7,file=outname,status='new' ) 
print* 

print*, '  Delta  t  =  ' 

print* 

n=0 

do  100  i=l, 10000 

read(4,903)  c(l) ,tim( 1) ,c(2) ,amp( 1) ,c(3) ,c(4) ,tim(2) , 

&  c( 5 ) , amp( 2) , c( 6) 

if  (tim( 1).  eq. 9. 999999)  then 

if  (amp( 1). eq. 9. 999999)  go  to  200 
endif 
n=n+l 

write(7,904)  tim( 1) ,amp( 1) 
if  (tim(2). eq.  9.  999999)  then 

if  (amp(2). eq.  9.  999999)  go  to  200 
endif 
n=n+l 

write(7,904)  tim(2) ,amp(2) 
if  (n. eq. 2)  go  to  300 
if  (n. eq. 102)  go  to  300 
if  (n.  eq.  202)  go  to  300 
if  (n. eq. 302)  go  to  300 
if  (n. eq. 402)  go  to  300 
if  (n. eq. 502)  go  to  300 
if  (n. eq. 602)  go  to  300 
if  (n. eq. 702)  go  to  300 
if  (n. eq. 802)  go  to  300 
if  (n. eq. 902)  go  to  300 

100  continue 
print* 

print*,'  There  are  more  points  remaining  in  HP  Vista  data  file' 
print*,'  increase  the  loop  size' 

200  write( 7 ,904)  9999.  ,9999. 
print* 

print*,'  Total  number  of  data  points  is' 
write(5,905)  n 
go  to  900 

300  dt=tim( 2) -tim( 1) 
print*, dt 
go  to  100 

900  close  (unit=4) 
close  (unit=7) 

901  format  (a20) 

902  format  (il) 

903  format  (3x,al,el5. 6,al,el5. 6,al,2x,al,el5. 6,al,el5. 6,al) 

904  format  (2x,el6. 8,5x,el6. 8) 

905  format  (lx,il0) 
end 
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APPENDIX  B.  DATA  TRANSFER  FROM  SOURCE  TO  VAX 

COMPUTER 


A.  DATA  SOURCE  TO  HP3565  COMPUTER 

Numerous  sources  of  data  were  used.  Analog  sources  included  function  generators 
and  FM  tape  recordings  of  accelerometer,  velocity,  and  pressure  transducers.  Another 
source  of  data  resulted  from  direct  measurements  of  vibrations  of  machinery  in  the  lab¬ 
oratory.  All  signals  were  fed  into  an  input  module  of  an  HP3565  computer  based  on  the 
HP9000  350  CPU  which  included  a  dynamic  signal  analyzer.  Hewlett-Packard  s 
HP- Vista  software  was  used  to  control  the  hardware.  The  input  module  was  calibrated 
and  ranged  for  the  applicable  voltages  of  the  input  signal.  The  digital  filters  were  set  to 
allow  the  frequencies  of  interest  to  pass.  The  primary  purpose  of  the  HP3565  computer 
system  was  to  digitize  and  filter  the  data. 

B.  HP3565  COMPUTER  TO  PC 

Once  the  data  was  digitized,  a  copy  of  the  raw  data  was  printed  out  for  use  in  veri¬ 
fying  the  data  upon  arrival  in  the  VAX  computer.  The  digitized  data  was  transferred 
from  the  HP3565  to  a  personal  computer  (PC)  as  a  printed  ASCII  file  output.  The  ac¬ 
tual  data  transfer  took  place  over  an  RS232  cable  connecting  the  laser  printer  port  on 
the  HP3565  to  the  PC  communications  port. 

C.  PC  TO  VAX 

At  the  PC  the  data  file  was  transferred  to  a  5.25"  floppy  disc.  This  data  was  down¬ 
loaded  to  the  VAX  computer  in  the  Mechanical  Engineering  CAD  CAE  Lab  via  the 
terminal  available  there. 

D.  EDITING  DATA  FILE  UPON  ARRIVAL  IN  THE  VAX 

The  data  file  which  arrived  in  the  VAX  included  a  banner  page  and  file  header  in¬ 
formation  which  had  to  be  manually  removed  by  editing.  In  addition,  there  were  form¬ 
feed  characters  throughout  the  file  which  had  to  be  removed.  The  data  in  the  VAX  was 
compared  with  the  data  printouts  made  by  the  HP3565.  On  the  average  there  were  3 
or  4  corrections  to  a  1024  data  point  file  due  to  the  high  data  transfer  rate  across  the 
RS232  cable  (9600  Baud).  A  slower  data  transfer  rate  could  have  been  used  to  improve 
the  accuracy,  but  the  transfer  time  would  have  taken  much  longer.  A  final  line  of 
characters  was  added  to  the  data  file  to  serve  as  an  end  of  file  indication.  A  computer 
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program,  CONVERT. FOR  (the  last  program  in  Appendix  A),  was  used  to  convert  the 
data  from  the  edited  IIP3565  format  into  the  format  used  in  subroutine 
DA  TAIN. INCLUDE.  The  data  format  conversion  program  also  provided  a  final  check 
to  ensure  that  all  data  points  had  been  included  by  counting  the  number  of  data  points. 
A  block  size  setting  of  1024  real  or  S192  real  on  the  HP3565  resulted  in  1024  or  8192 
data  points,  respectively. 
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